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Chapter  1 


Introduction 


Finite  element  methods  use  polynomials  to  approximate  the  solution  to  boundary 
value  problems  over  some  subdivision  (mesh)  of  the  domain.  The  traditional  method, 
called  the  h  version,  involves  keeping  the  degree  of  the  approximating  polynomial 
fixed,  usually  at  a  low  level,  and  achieves  accuracy  by  properly  refining  the  mesh. 
A  more  recent  method  is  the  p  version,  where  the  mesh  is  fixed  and  accuracy  is 
achieved  by  increasing  the  degree  of  the  polynomial  approximation.  The  hp  version 
is  a  combination  of  these  two  approaches. 

Currently,  there  are  several  commercial  finite  element  codes  available  or  under 
development  including  RASNA,  PHLEX,  STRESSCHECK,  and  IBM-POLYFEM  that 
offer  these  new  p  and  hp  methods.  Moreover,  codes  well  established  in  industry,  such 
as  MSC/NASTRAN,  are  being  equipped  with  p/hp  technology. 

Codes  with  hp  capabilities  allow  the  user  to  selectively  employ  a  combination  of  h- 
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refinement  and  p-refinement  to  achieve  accuracy.  For  example,  a  highly  refined  mesh 
with  low  p  may  be  used  near  corners,  while  large  elements  with  high  p  may  be  used 
in  the  interior.  The  selection  of  h  and  p  may  also  be  done  adaptively. 

The  increased  flexibility  afforded  by  such  hp  codes  puts  increased  demands  on  the 
underlying  variational  method  and  elements  being  used.  The  code  must  now  be  robust 
with  respect  to  both  h-refinement  and  p-refinement.  Various  quantities  of  engineering 
interest,  in  particular  the  displacements  and  stresses,  should  be  calculated  robustly 
and  optimally  for  the  type  of  refinement  chosen.  The  correct  design  of  elements  and 
underlying  method  is  therefore  crucial. 

For  many  materials  and  geometries,  this  level  of  robustness  is  easily  realized  by 
using  the  standard  finite  element  method  (SFEM).  However,  in  the  case  of  isotropic 
linear  elasticity,  if  the  material  is  almost  incompressible  (Poisson  ratio  u  close  to  ^), 
the  SFEM  can  lead  to  dilational  (or  Poisson  ratio)  locking  (see  e.g.  [1,  6,  13,  25]). 
Locking  is  when  the  rate  of  convergence  decreases  as  In  severe  cases  the  method 

may  not  converge  to  the  correct  solution  at  all.  This  is  an  important  exception  since 
almost  incompressible  materials  are  commonly  used  in  industry.  Natural  rubber  is 
almost  incompressible  and  materials  that  undergo  plastic  deformations  may  also  be 
considered  almost  incompressible  [23].  Certain  biological  tissues  are  also  considered 
nearly  incompressible.  Also,  the  limiting  case  of  z/  =  |  reduces  to  Stokes  problem, 
which  is  obviously  of  great  importance  in  fluid  flow  problems. 

In  the  context  of  isotropic  linear  elasticity  under  conditions  of  plane  strain,  it 
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is  well  known  from  computational  experience  that  locking  occurs  in  various  finite 
element  schemes.  For  example,  with  v  close  to  0.5,  the  h  version  using  piecewise 
linear  polynomials  on  triangular  meshes  produces  a  very  poor  approximation  of  the 
displacements.  However,  in  1985,  Scott  and  Vogelius  (see  [18]  and  [19])  showed  that 
the  h  version  is  locking  free  for  both  displacements  and  pressures  ^  when  polynomials 
of  degree  A:  >  4  are  used  on  certain  triangular  meshes.  ^  They  showed  that  locking 
can  also  be  avoided  for  A;  <  4  if  very  specialized  meshes  can  be  constructed.  In  1992, 
Babuska  and  Suri  [1]  characterized  displacement  locking  for  A;  <  4,  which  was  shown 
to  depend  not  only  on  k  but  also  on  the  mesh.  They  also  showed  that  for  two  types 
of  rectangular  elements,  locking  cannot  be  avoided  for  any  k. 

The  p  version  has  long  been  known  to  be  locking-free  for  displacements.  In  1983, 
Vogelius  [26]  showed  that  for  computing  displacements  on  triangular  meshes,  the  p 
version  converges  uniformly  in  i/,  with  optimal  rate  up  to  an  arbitrary  e  >  0.  Babuska 
and  Suri  [1]  removed  the  dependence  on  e  in  the  locking-free  rate  of  convergence  for  the 
p  version.  They  also  showed  this  applies  to  parallelogram  elements  as  well.  However, 
the  references  cited  here  do  not  guarantee  the  pressures  will  be  locking  free  for  the  p 

version  (though  results  exist  for  the  h  version). 

^  We  use  the  term  “pressures”  in  place  of  stresses  here  because  the  quantity  that  suffers  most  due 
to  locking  is  the  sum  of  the  normal  stresses,  which  reduces  to  a  multiple  of  the  pressure  in  Stokes 
problem,  i.e.  in  the  limit  v  —  (See  equation(2.6).)  We  use  the  two  terms  interchangeably  in  the 
sequel. 

^ Since  p  will  be  used  to  denote  pressure,  we  use  k  to  denote  polynomial  degree. 
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Mixed  methods  have  also  been  used  to  avoid  locking.  In  these  methods,  an  aux¬ 
iliary  variable  related  to  the  stresses  is  introduced,  and  the  solution  is  found  as  a 
saddle  point,  rather  than  the  minimum  of  the  energy  functional.  Several  mixed  p- 
type  elements  have  been  studied  for  the  Stokes  problem  in  the  context  of  the  spectral 
element  method [4].  If  chosen  appropriately  they  converge  well  in  both  displacements 
and  pressures.  This  is  also  true  for  h  version  mixed  methods.  The  difficulty  lies  in 
choosing  the  mixed  method  spaces  correctly.  As  shown  in  [10],  if  the  spaces  are  not 
properly  balanced,  mixed  methods  may  converge  at  less  than  optimal  rates,  or  they 
may  may  not  converge  at  all.  It  is  well  known  [6]  that  mixed  method  spaces  must 
satisfy  an  inf-sup  (Babuska  -  Brezzi)  condition  to  be  stable  and  that  choosing  such 
spaces  can  be  difficult.  More  specifically,  identifying  spaces  and  proving  that  they  are 
both  stable  and  nearly  optimal  can  be  difficult. 

The  problem  of  hp  mixed  methods  has  been  addressed  by  Stenberg  and  Suri[20]. 
They  identify  sufficient  conditions  for  selecting  mixed  method  spaces  on  parallelogram 
elements  that  produce  nearly  optimal  and  stable  methods.  Using  these  conditions,  they 
identify  several  mixed  methods.  They  show  that  these  mixed  method  elements  are 
optimal  in  both  displacements  and  pressures  for  h-refinement,  and  that  jo-refinement 
produces  at  worst  0{k^)  locking  in  displacements  and  locking  in  pressures. 

However,  this  study  did  not  include  any  niimerical  verification  of  these  methods. 

In  the  p  version,  since  the  element  size  remains  fixed,  curved  elements  are  essential 
for  practical  problems.  In  [22],  some  experiments  with  curved  elements  showed  that 
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both  the  h  and  p  versions  of  SFEM  show  increased  locking  in  the  displacements  when 
curved  elements  are  used.  For  mixed  h  version  methods,  the  Q2  —  Vo  element  (see 
[6])  is  known  to  be  stable  for  curved  elements,  as  are  some  other  elements. 

Our  goal  in  this  dissertation  is  to  identify  finite  element  methods  that  compute  both 
displacements  and  stresses  with  uniform,  near  optimal  accuracy  for  any  u  G  [0,  |), 
when  either  h-  or  p-refinement  is  used.  In  addition,  the  methods  should  provide 
nearly  optimal  performance  on  domains  of  practical  interest,  which  include  curvilinear 
domains. 

Several  definitions  and  descriptions  are  needed  throughout  this  dissertation;  these 
are  given  in  Chapter  2.  We  describe  the  differential  and  variational  forms  of  linear 
elasticity  and  the  associated  finite  element  methods.  We  then  define  locking  and  related 
issues. 

In  Chapter  3  we  do  a  computational  study  of  finite  element  methods  on  paral¬ 
lelogram  elements.  We  summarize  the  available  theory  mentioned  above  so  we  can 
compare  this  with  our  computational  results.  We  examine  the  mixed  methods  from 
[20]  to  see  how  they  perform  numerically.  In  [20] ,  it  was  stated  that  the  inf  —  sup 
constant  for  these  methods  decayed  at  a  rate  no  worse  than  k~^  —  we  validate  this 
numerically.  The  conditions  used  in  [20]  to  define  mixed  methods  indicate  that  a  tra¬ 
ditional  spectral  element,  ,  is  not  an  optimal  choice.  This  is  because  the 

polynomial  space  for  pressures  should  be  one  degree  less  than  the  displacements,  and 
not  two.  We  computationally  compare  this  element  with  mixed  methods  from  [20] 


6 


and  show  that  the  element  is  indeed  less  than  optimal.  We  also  show  that 

for  a  model  problem  these  methods  work  well  for  computing  both  displacements  and 
pressure.  The  STEM  is  also  compared  with  these  mixed  methods  . 

We  turn  our  attention  to  straight  sided  triangular  elements  in  Chapter  4.  The 
results  in  [26]  imply  that  the  divergence-stability  constant  for  the  p  version  SEEM 
on  triangular  elements  is  k~°‘  for  some  unknown  o;  >  0.  Since  triangular  elements 
are  used  extensively  in  several  hp  commercial  codes,  we  investigate  their  performance 
computationally  as  z/  — >  0.5.  We  find  that  the  stresses  can  be  extracted  fairly  accur¬ 
ately,  which  suggests  that  a  in  the  estimate  from  [26]  is  small.  However,  when  we 
evaluate  curved  triangular  elements,  we  find  that  stress  extraction  is  less  satisfactory. 
We  therefore  seek  mixed  methods  for  which  a  better  characterization  of  the  inf  —  sup 
constant  can  be  made.  These  mixed  triangular  elements  are  formulated  to  satisfy  the 
conditions  outlined  for  parallelograms  in  [20].  One  of  these  elements  was  analyzed  in 
[17],  where  a  lower  bound  of  k~^  was  obtained  on  the  inf —  sup  constant.  We  show 
that  this  bound  of  k~^  is  quite  pessimistic  in  the  range  of  k  used  in  practice  by  (a) 
computing  the  inf  —  sup  constant  numerically  and  (b)  testing  the  rate  of  convergence 
computationally  for  some  test  problems.  Let  us  mention  that  the  elements  we  con¬ 
sider  use  discontinuous  pressures.  In  the  case  of  a  continuous  pressure  space,  Boillat 
[5]  showed  that  the  divergence-stability  constant  behaves  like  k~'^  ^\/ln  for  the  p 
version  based  on  Taylor-Hood  elements. 

Since  curved  elements  are  essential  for  p  version  meshes,  in  Chapter  5  we  invest- 
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igate  locking  on  curved  elements.  We  prove  that  on  a  single  curved  element,  the 
displacement  error  behaves  like  (k  —  where  a  depends  on  the  mapping  from  a 
reference  element  to  the  curvilinear  element.  Next,  we  consider  the  mixed  methods 
from  [20]  on  curved  elements  and  show  how  the  stability  constant  can  be  estimated 
for  curvilinear  meshes.  We  use  this  result  to  compute  inf  —  sup  constants  for  several 
maps  and  show  that  they  are  stable  for  practical  elements.  We  also  formulate  a  new 
curvilinear  quadrilateral  element  and  show  that  the  inf-sup  constant  for  this  element 
is  bounded  below  by  Ck~2.  This  is  the  same  bound  as  for  parallelogram  elements. 

When  the  pressure  is  of  interest,  we  show  that  the  p  version  SFEM  shows  significant 
deterioration  when  curved  elements  are  used.  We  show  experimentally  that  this  effect 
is  worse  when  elements  are  curved  at  the  boundary  as  opposed  to  having  curved 
interior  edges.  We  also  perform  experiments  to  show  that  the  mixed  methods  from 
Chapters  3  and  4  are  robust  for  a  model  problem. 

As  a  final  test,  in  Section  5.4,  we  use  a  benchmark  problem  from  [25]  to  invest¬ 
igate  the  performance  of  SFEM  and  a  mixed  method  for  computing  point  values  of 
stresses.  In  [25],  Szabo,  Babuska  and  Chayapathy  showed  that  SFEM  is  ineffective  for 
computing  point  values  of  the  sum  of  normal  stress  and  proposed  a  post-processing 
method  to  overcome  this  shortcoming.  We  demonstrate  that  curvilinear  mixed  method 
elements  can  be  used  to  accurately  extract  these  point  values. 

Recently,  Noel  and  Szabo  [16]  proposed  a  spatial  formulation  for  geometrically 
nonlinear  elasticity.  They  showed  that  their  SFEM  works  well  on  many  problems  when 
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u  is  not  close  to  0.5.  In  Chapter  6,  we  reformulate  this  problem  with  a  mixed  method 
and  investigate  locking  of  the  standard  and  the  mixed  method.  We  find  that  while  the 
standard  method  produces  erratic  results  when  computing  point  values  of  stresses, 
the  mixed  method  eliminates  the  oscillations.  This  result  mirrors  that  of  Section  5.4 
for  a  linear  problem.  Our  results  therefore  suggest  that  the  mixed  methods  discussed 
here  can  be  applied  to  eliminate  locking  in  such  non-linear  problems  as  well. 

The  overall  conclusion  of  this  thesis  is  that  the  mixed  methods  investigated  here 
are  excellent  candidates  for  hp  implementation  when  a  robust  code  is  desired. 

Notation 

In  this  dissertation,  we  shall  denote  the  n-dimensional  Euclidean  space  by  R",  with 
X  =  (a;i,a;2)  or  x  =  {x,y)  and  |x|  =  {x^  +  We  also  denote  the  set  of  non¬ 

negative  integers  by  N.  If  Q  is  any  one-  or  two-dimensional  set,  Q  denotes  its  closure. 
By  we  denote  a  domain  in  R^  with  boundary  dQ.  By  (fi)  we  denote  the  Sobolev 
space  of  functions  on  0  with  square  integrable  generalized  derivatives  of  order  <  r 
(r  >  0)  furnished  with  the  norm 

0<|a|<r 

where  a  =  (q:i,Q!2)  7  >  0,  f  =  1,2  integers,  |a|  =  ai  -1-  02,  and 

f^^TL  z=z  - 

dx'l^dx^^' 
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As  usual,  (0)  =  (fl).  For  functions  u,v  ^  (0) ,  {u,v)  will  denote  the  usual 

inner  product.  We  denote 

Ll  (0)  ^  !^U  e  (n)  :  J^udx  =  oY 

Also, 

Hq  (0)  =  |u  G  (O)  :  u  =  0  on  j- , 
and  we  will  further  use  the  semi-norm  notation 

|a|=r 

Throughout  this  dissertation,  the  letter  C  will  be  used  to  denote  generic  positive 
constants,  possibly  not  the  same  in  each  occurrence.  Finally,  K  will  refer  to  an 
element  in  a  finite  element  mesh  Ch  on  a  domain  ft. 


Chapter  2 
Preliminaries 

2.1  The  Standard  and  Mixed  Methods  for  Linear 
Elasticity 

The  differential  form  of  the  equations  of  isotropic  linear  elasticity  consists  of  three 
relationships.  They  include  the  equilibrium  equations 

(2.1a)  -(Tijj  =  fi  on  P, 

the  constitutive  relation  for  isotropic  plane  strain 

(2.16)  aij  =  Adivu^ij  +  2fieij{u), 

and  the  deformation  (strain)  tensor 

10 


(2.1c) 
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which  relates  the  displacement  u  =  {ui,U2)  to  the  strain.  We  let  u  =  0  on  To  and 
cTijUj  =  gi  on  Fi,  where  =  FoUFi  and  FonFi  =  0.  Let  0  <  z/  <  0.5  be  the  Poisson 
ratio  and  E  the  modulus  of  elasticity,  and  define  the  Lame  constants 

^  ^  ^  Eu  E 

(l  +  z/)(l-2i/)’  ^“2(1+1/)' 

Then  the  variational  form  of  our  problem  is:  Find  u  G  V  such  that  for  all  v  G  V, 
(2.3)  (e  (u) ,  e  (v))  +  A  (div  u,  div  v)  =  F  (v) , 

where 

F  (v)  =  j J  f  -vdx  +  g-vds 

and  V  =  [H^  For  the  case  Fq  =  0,  we  assume  that  F  {TZ)  =  0  for  any  rigid  body 

motion  Tl.  This  ensures  that  (2.3)  has  a  unique  solution  (modulo  rigid  body  motions 
—  these  are  assumed  eliminated  in  computations  by  suitable  constraints). 

Given  a  sequence  of  finite  element  subspaces  {V  n}  N  C  V,  we  define  the  stand¬ 
ard  finite  element  approximation  to  (2.3)  as:  Find  G  Vw  such  that  for  all  v  G  V^r, 


(2.4)  2fi  (£  (uiv) ,  £  (v))  +  A  (div  u^v,  div  v)  =  F  (v) . 

We  identify  the  parameter  N  with  the  dimension  of  the  subspace  Vat-  Associated  with 
this  standard  method,  we  define  the  energy  norm 


(2.5) 


u|l£:,^=  (2/xll£(u)  l|2-|-Alldiv  (u)ll^)"  . 
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Next,  we  describe  the  mixed  formulation.  We  define  the  new  independent  unknown 

(2.6)  Adivu, 

which  is  a  multiple  of  the  sum  of  the  normal  stresses,  i.e.  p  =  +  ^22)  = 

2  (I:^)  +  ^y)  ^  Ihis  corresponds  to  the  pressure  in  the  limiting 

Stokes  equations.)  Then  (2.3)  is  written  in  the  Herrmann  variational  form  [9]:  Find 
(u,p)  G  V  xW  such  that  for  all  (v,  q)  G  V  x  W, 

(2.7)  2p  (e  (u) ,  e  (v))  -  (p,  div  v)  =  F  (v) 

(2.8)  (div  u,  q)  +  ^  (p,  g)  =  0 

where  W  =  (0). 

For  the  mixed  finite  element  method,  we  assume  we  are  given  a  sequence  of  finite 
element  subspaces  (Vjv  x  Wn)  C  V  x  VF  and  find  {un,Pn)  £  (V^r  x  kFv)  such  that 

(2.9)  2p  (e  (uat)  ,  e  (v))  -  {pN,  div  v)  =  F  (v) 

(2.10)  (div  uat,  9)  +  j  {pn,  q)  =  0 

for  all  (v,  q)  G  Vjv  X  Wn-  The  correct  analog  of  the  “energy  norm”  for  mixed  methods 
should  be 

(2.11)  \Me,u^  (2/r||e(u)||^  + A~^||p||^)"  . 

As  is  well  known  [6],  the  accuracy  of  this  method  will  not  only  depend  upon 
how  well  VjVjHiv  approximate  V,  VF  respectively,  but  also  the  stability  of  the  pair 


^Note  that  an  =  a^  and  (T22  =  ay. 
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(VnjWn)  (i.e.,  the  inf-sup  or  Babuska-Brezzi  condition  satisfied  by  them).  This  is 
made  precise  in  the  following. 


Theorem  2.1  [6]  Suppose  (Vjv  x  IT/v)  C  (V  x  W)  and  let  (u,p)  satisfy  (2.7)  and 
(2.8)  and  let  (uniPn)  satisfy  (2.9)  and  (2.10).  If  there  exists  /?  >  0  such  that 


(2.12) 


inf  sup 

qeWtf  veVjv 


l(gidivv) 

lkllo||v||i 


>f3 


then 


u;v||i  +  lb  -  P^llo  <  C  Qnf^  l|u  -  v||i  -f  \\p 


2.2  Locking 

Before  we  give  the  formal  definitions  relative  to  locking,  we  explain  the  practical 
meaning  by  way  of  an  example.  In  (2.4)  we  see  that  if  A  — )■  oo,  the  term  (divu;v^,  divv) 
must  approach  zero  for  the  equation  to  have  meaning.  This  forces  divuiv  — >■  0.  When 
the  finite  element  space  is  not  rich  enough  to  simultaneously  satisfy  this  constraint 
and  approximate  u  well,  we  can  expect  poor  performance.  In  linear  elasticity,  A  — )•  oo 
has  physical  meaning;  it  is  the  case  of  nearly  incompressible  materials. 

We  illustrate  what  can  happen  as  A  — >•  oo  in  Figure  2.1,  which  duplicates  Figure  4.2. 
This  shows  the  results  of  using  the  h  version  SFEM  to  solve  (2.4)  on  a  uniform 
triangular  mesh.  Recall  that  A  is  related  to  u  by  (2.3),  and  A  — )■  oo  is  equivalent  to 
z/  — We  see  that  when  k  =  2  the  rate  of  convergence  in  energy  norm  changes  from 
0{h^)  when  u  —  0.3  to  0(h)  when  u  =  0.4999.  In  contrast,  when  A;  =  4,  the  energy 
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Degrees  of  Freedom  —  h-refinement  Degrees  of  Freedom  —  h-refinement 


Figure  2,1:  ST,  /j-refinement,  Smooth  Solution. 

norm  error  curves  for  u  =  0.3  and  0.4999  are  nearly  identical.  This  deterioration  in 
the  energy  norm  rate  of  convergence  as  i/  |  is  known  as  locking.  Notice  that  locking 
also  occurs  when  A;  =  2  for  the  SNS  (sum  of  normal  stresses)  relative  error.  When 
u  =  0.4999  the  error  curve  actually  increases  as  h  decreases.  Another  important  effect 
is  observed  when  A:  =  4  for  the  SNS  error.  The  curves  for  v  =  0.3  and  0.4999  have  the 
same  slope,  but  the  u  =  0.4999  curve  is  shifted  up  relative  to  v  —  0.3  curve.  This  shift 
at  a  given  level  of  refinement  is  referred  to  as  the  locking  ratio.  A  desirable  feature  of 
a  finite  element  method  is  to  be  free  of  locking  and  to  have  a  locking  ratio  near  unity 
over  the  range  of  discretizations  of  interest. 

We  now  make  these  concepts  more  precise.  We  will  need  the  following  ingredients 
for  onr  definitions  [22]. 

1.  Solution  space  Hj,.  We  assume  that  our  exact  solutions  s,,  lie  in  some  sets 
H,,,  which  characterize  the  smoothness  of  s^.  For  solutions  Sj,  of  comparable 
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smoothness  we  expect  comparable  results,  so  we  will  define  locking  with  respect 
to  the  family  {H^}.  We  assume  that  ||  •  ||h^  is  a  norm  on  Hj,  and  define  for  any 
B  >  0,  the  set 

Hf  =  {s,  e  H,  :  IIs.IIh.  <  B}. 

2.  Error  functional  E^,.  Next,  we  assume  we  are  interested  in  a  certain  functional 
Ej,  :  Sy  ->  R.  We  will  be  interested  in  the  norm  of  the  displacement  and 
the  norm  of  the  sum  of  normal  stresses.  We  should  point  out  that  locking 
is  closely  coupled  to  the  error  functional.  A  method  may  be  locking  free  in  one 
error  measure  but  show  significant  locking  in  another. 

3.  Extension  procedure  S.  We  assume  there  is  a  sequence  of  finite  element  spaces 
^  =  {Vjv}-  Then  S  gives  an  extension  procedure,  i.e.  a  rule  on  how  to  increase 
the  dimension  N  (e.g.  by  the  h  version,  p  version,  etc.).  We  assume  here  that 
N  ranges  over  a  set  Af.  As  we  shall  see,  locking  is  strongly  dependent  on  8. 

4.  Parameter  set  S.  Our  parameter  here  is  u,  the  Poisson  ratio,  which  we  assume 

varies  over  the  set  5  =  {i/  :  0  <  <  0.5}.  (In  computations,  we  will  often  take 

S  =  {0.3,0.4999}. 

We  first  define  locking  ratio. 


Definition  2.1  The  function  L{u,  N),  called  the  locking  ratio  for  (u,  N)  G  N  x  Af,  is 
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given  by 

l{y  N)  =  -  O 

inf^e5„  sup^^gHB  E^{w^  -  w^)  ’ 
where  Sa  =  S  Q  (—00,  a]  for  some  a  <  0.5  such  that  Sa  ^  0- 

Obviously,  L  also  depends  on  and  Sa-  It  compares  the  performance  of 

the  method  at  Poisson  ratio  u  to  the  best  possible  performance  for  reasonable  values 
of  u,  characterized  hy  v  <  a. 

Let  us  now  present  the  definition  of  locking  ([2],  [1]),  for  the  case  that  a  >  0. 

Definition  2.2  The  extension  procedures  is  locking  free  for  u  G  [0,0.5),  with  respect 
to  the  solution  sets  and  error  measures  E^,  if  and  only  if 

limsup  I  sup  L(i/,  N)\  =  C  <  00. 

N^oo  \i/6[0,0.5)  / 

S  shows  locking  of  order  f(N)  if  and  only  if 
0  <  lim  sup  j  sup 

N-^00  \  V 

where  f{N)  00  as  N  00. 

The  denominator  in  (2.13)  can  be  replaced  by  Eo(N),  the  asymptotic  rate  of  best 
approximation  in  of  functions  in  by  functions  in  Vn, 

Eo{N)=  sup  inf  l|u;-u||i,n, 

which  represents  the  smallest  error  which  could  be  achieved  using  Vn  to  approximate 
the  most  unfavorable  w  in  (i/^)®.  For  example,  for  the  h  version,  Eq{N)  =  0{h^)  = 


L{u,  N) 


f{N) 


C  <  00 
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0{N  t)  where  =  min{A:  —  l,p},  while  for  the  p  version,  Fo[N)  =  0{p  = 

0{N-^). 

For  the  standard  method,  we  are  interested  in  two  error  functionals: 

1.  The  norm  of  the  displacement  error  (which  is  equivalent  to  the  energy  norm 
(2.5)). 

2.  The  norm  of  the  error  in  the  sum  of  normal  stresses,  SNS  (which  is  computed 
by  differentiating  the  displacement). 


For  the  mixed  method,  we  use  the  same  error  measures,  but  we  compute  the  energy 
norm  using  (2.11)  and  we  compute  SNS  directly,  since  it  is  a  scalar  multiple  of  the 
independent  variable  p. 


Chapter  3 


Parallelogram  Elements 


We  describe  the  SEEM  and  several  mixed  methods  [20]  on  parallelogram  elements.  We 
show  experimentally  that  for  these  mixed  methods  the  inf  —  sup  constant  does  indeed 
decay  like  k~^.  We  then  simimarize  the  rates  of  convergence  for  these  methods.  For 
a  model  problem,  we  show  that  these  methods  are  locking  free  for  computing  both 
displacements  and  pressure.  We  also  compare  the  SFEM  with  these  mixed  methods. 

3.1  Approximation  Spaces 

3.1.1  General  conditions  for  stability  and  approximability 

We  begin  with  the  notation  and  definitions  needed  to  define  the  finite  element  spaces 
we  will  use.  We  do  this  here  not  only  for  parallelogram  elements,  but  for  the  other 
elements  considered  in  this  thesis  as  well. 
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Assume  that  a  sequence  of  meshes  {Ch}  consisting  of  (curvilinear)  quadrilaterals 
and/or  triangles  is  defined  on  O,  which  is  regular  in  the  following  sense  [7]: 

i)  The  intersection  of  two  elements  in  Ch  consists  of  one  point  (a  vertex),  one  entire 
common  edge  or  the  empty  set, 

ii)  There  exists  a  constant  a  such  that 

V/T  eCh,—<  a. 

PK 

Here,  hf[  is  the  diameter  of  K  and  px  is  the  diameter  of  the  largest  inscribed  ball  in  K. 
Ch  does  not  have  to  be  quasi-uniform  (see  Remark  3.3  below).  For  K  €  Ch,  we  denote 
to  be  the  smooth  invertible  mapping  from  the  reference  element  K  onto  K,  where 
K  =  [—1, 1]^  when  K  is  a  quadrilateral  and  K  =  {{x,y)  :  —I  <  x  <  1,-1  <  y  <  —a:} 
when  K  is  a  triangle.  (Further  restrictions  on  J-x  will  be  imposed  in  the  chapter  on 
curvilinear  elements.) 

Let  'VkiK)  and  Wk{K)  be,  respectively,  the  spaces  for  u  and  p  on  K,  consisting 
of  polynomials  of  degree  related  to  k.  Define  for  any  K  £  Ch, 

(3.1)  Yh{K)  =  {v  =  vo:f-i  :  vG  Vfc(A0}, 

Wk{K)  =  {p^poJ^J,^  ■.  p£Wk{k)]. 

Then  with  N  =  N{h,  k),  the  finite  element  spaces  are  defined  as 


(3.2) 


Vjv  =  {vG  V  :  V|kG  Vfc(/T)  VA^GC,}, 

Wn  =  {pew  :p^xeWk{K)\/K£Ch}, 
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where  V  and  W  are  the  solution  spaces  for  (2.7)  and  (2.8).  We  will  denote  by  Vk{K) 
the  set  of  polynomials  of  total  degree  k  on  while  Qk{K)  will  be  the  set  of  poly¬ 
nomials  of  degree  k  in  each  variable.  The  serendipity  or  trunk  space  Q'y.{K)  will  be 
span  {Pkik), 

As  noted  earlier,  the  mixed  method  spaces  we  consider  are  based  on  [20].  In 
that  paper,  six  conditions  are  given  for  selecting  mixed  method  spaces  Vi;,Wfc  on 
parallelogram  elements.^  These  conditions  motivate  the  choice  of  the  mixed  finite 
element  spaces,  and  are  central  to  arguments  developed  in  later  chapters.  Since  we  are 
interested  in  extending  the  mixed  methods  of  [20]  to  more  general  meshes,  (including 
straight  sided  triangular  and  curvilinear  triangular  and  quadrilateral  meshes),  we  will 
list  a  more  general  form  of  two  of  these  conditions. 

The  first  condition  in  [20]  is  that  identical  subspaces  are  used  for  both  displacement 
components,  i.e., 

(Ai)  \,(k)  =  [Vt(k)]\ 

Let  V^{K)  =  Vk{K)  n  Hq{K)  denote  the  set  of  internal  shape  functions  used  for 
all  components  of  the  displacement,  so  that  V^(.^)  =  Yk{K)  fl  [Hq(K)]^  =  [V^{K)Y. 
We  note  that  there  exists  a  space  Xk{K)  such  that 

V^{k)  =  {u|u  =  w  G  Xk{k)}, 

where  is  the  lowest  degree  “bubble  function”  on  K,  the  reference  element  (square 


^The  conditions  from  [20]  are  given  for  both  and  We  restrict  our  presentation  to 
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or  triangle).  The  second  condition  is  that  the  space  for  the  pressure  satisfies  the 
following. 

(A2)  Vq  G  [Xkik)]^  V  q  €  Wk{k). 

Due  to  the  definitions  of  the  spaces  and  Xk{K),  we  can  define  a  weighted 
projection  Tk  :  Hoi^)  by 

(3.3)  {v~rkv,w)j^  =  0  \/weXkik). 

The  third  condition  is  that  the  projection  operator  Tit  satisfies 

(A3)  \\Tku\\-^^  fr  <  for  some  7  >  0. 

Note  that  in  [20],  instead  of  a  general  7,  condition  (A3)  is  stated  with  7  =  |,  since 
for  parallelogram  elements  considered  there,  this  can  be  proven. 

Conditions  (A1)-(A3)  give  a  local  stability  condition  for  any  element  K  for  which 
J-ji  is  linear.  The  following  lemma  follows  from  Lemma  5.1  of  [20]. 

Lemma  3.1  Let  the  spaces  Yk{k),Wk{k)  satisfy  conditions  (Al)-(A3).  Let  K  = 
J-K{k),  where  Tk  is  an  affine  mapping.  Then  for  every  q*  G  Wk{K)  fl  Lq{K)  there 
exists  V*  G  Yk{K)  n  such  that 

(3.4)  (divv*,g*)A'  >  Cik~'^\\q’^\\lK  and  \v*\i,k  <  C2\\q*\\o,K-, 

where  (71,(72  are  positive  constants  independent  of  K,h,k  and  q*. 

Note  that  7  in  (3.4)  may  be  different  over  different  K.  In  this  case,  the  value  chosen 
in  (3.4)  is  the  maximum. 
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Remark  3.1  Let  us  remark  that  if  (Al)  and  (A2)  hold,  then  we  can  always  define 
Tk  by  (3.3).  For  V^{K)  ^  0^  we  will  always  have 

for  some  C(k)  <  oo.  Then  Lemma  3.1  holds  with  k~'^  replaced  by  (C(k))~^ .  Hence, 
(Al),  (A2)  are  sufficient  to  ensure  a  unique  solution  to  (2.9)-(2.10). 

As  we  shall  see  in  later  chapters,  the  best  theoretical  estimate  available  for  7 
in  (AS)  is  pessimistic  for  triangular  elements,  and  unavailable  for  most  curvilinear 
elements.  In  order  to  investigate  such  elements,  we  will  estimate  the  local  inf-sup  con¬ 
dition  (3.4)  computationally.  Once  this  local  inf-sup  constant  has  been  estimated  for 
all  elements  in  a  given  mesh,  the  remaining  theory  developed  in  [20]  is  still  applicable. 

We  therefore  formulate  the  following  alternative  condition  to  (A1)-(A3). 

(LS)  For  every  K  €  Ch,  given  q*  G  Wk{K)  fl  Lq{K),  there  exists  v*  G  Yk(K)  fl 
such  that  (3.4)  holds. 

Let  q  G  Wn  be  arbitrary  and  write  q  =  q  +  q*  with  q  being  the  L^  projection  of  q 
onto  the  space  of  piecewise  constants: 

WN  =  {q^WN:  q\K  €  Vo{K)  MK  G  Ch}. 

To  have  a  global  condition  we  need  the  following. 

(A4)  The  pair  (Vjv,  Wn)  is  stable, 
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where  Vjv  =  {v  :  V|/f  =  v  o  for  some  v  G  'E2{K)  WK  G  Ch},  where  Fi2{K)  = 
[Q2{K)^  if  K  is  the  reference  square  and  'Ei2{K)  =  \V2{K)Y  if  K  is  the  reference 
triangle.  Then  we  have  the  following,  see  [20,  Theorem  5.1]. 

Theorem  3.1  Let  the  spaces  YmiWn  satisfy  conditions  (A1)-(A4)  or  (LS),(A4). 
Then 

(3.5)  sup  >  Ck-'  ||,||„  V,  e  ITk, 

vev^,\{o}  l|vlli 

where  the  constant  C  is  independent  of  h,k  and  q. 

The  proof  is  similar  to  that  of  Theorem  5.2. 

Finally,  we  will  take  the  spaces  YkiK)  and  Wk{K)  to  contain  polynomials  of  degree 
k  and  fc  —  1  —  n,  respectively.  Note  that  in  order  to  get  the  optimal  convergence  rate 
in  both  h  and  fc,  we  require  n  =  0.  More  precisely,  we  need  the  following  conditions. 

(A5)  {Qk\K)Y  C  Yk{K)  if  K  is  the  reference  square, 

[Pk{K)]'^  C  Yk{K)  if  K  is  the  reference  triangle. 

(A6)  Vk-i-n{K)cWk{K),  n>0. 

We  will  assume  that  k  >2,  i.e.  we  exclude  the  case  with  piecewise  constant  pressure. 

3.1.2  Choice  of  Parallelogram  Spaces 

We  now  define  our  methods  by  fixing  the  choices  of  Vk{K)  and  Wk{K),  when  is  a 
parallelogram.  We  define  the  polynomial  spaces 
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Qk{k\Ql{k)cQu{k)hy 
Qk{k)  =  Qk{k)f\Vk+2{k\^.n^, 

Ql{k)  =  Qk{k)  augmented  by  all  bubble  functions  of  degree  A:  +  1,  i.e.  Qk{k) 
together  with  functions  in  Qk+i{k)  that  vanish  on  dk. 

We  then  define  the  following  methods  [20]: 

MQl:  Vk{k)  =  [Qk{k)]\  Wk{k)  =  Vk-i{k). 

MQ2:  Vk{k)  =  [Qk{k)]\  Wkik)  =  Vk-i{k). 

MQ3:  v,{k)  =  [Q’,^2(k)]\  Wkik)  -  Vk-lik). 

MQ4:  Vk{k)  =  [Qlik)]\  Wkik)  =  Qk-lik). 

MQ5:  Vkik)  =  [Qkik)]\  Wkik)  =  Qk-2ik). 

MQ6:  Vkik)  =  [Qkik)]\  Wkik)  =  Qk-2ik)  U  Vk-iik). 

In  addition,  we  define 

SQ:  U(/0  =  [%(*■)]". 

The  mixed  method  spaces  MQl  and  MQ5  are  illustrated  for  A;  =  5  in  Figure  3.1, 
where  spanning  sets  for  W^ik)  and  one  component  of  V^ik)  are  shown. 

We  note  that  given  Wkik)  =  Vk-iik)  the  minimal  displacement  space  satisfying 
(A2)  is  ,  the  displacement  space  for  MQl.  Comparing  MQ2  or  MQ3  with 

MQl,  we  see  that  the  only  difference  is  that  this  minimal  space  has  been  replaced  by 
a  larger  space.  This  increases  the  number  of  degrees  of  freedom  without  increasing 
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MQl  MQ5 


0  a;  5  0  a;  5 


Figure  3.1:  Polynomial  spaces  for  MQl  and  MQ5  when  k  =  5. 

the  asymptotic  approximation  rate.  (For  instance,  Qk{K)  has  dimension  {k  +  1)^ 
compared  to  +  7/s)/2  for  Qk{K).) 

On  the  other  hand,  MQ6  can  be  regarded  as  an  optimized  version  of  MQ5, 
where  Vk{K)  has  been  kept  fixed  and  the  maximal  space  chosen  for  Wk{K).  This 
only  involves  adding  two  extra  degrees  of  freedom  per  element,  but,  as  Theorem 
3.4  below  shows,  leads  to  an  extra  order  in  the  asymptotic  rate  of  h  convergence. 
Similarly,  in  MQ4,  the  minimal  space  has  been  selected  for  Vk{K)  with  the  choice 
Wk{K)  =  Qk-i{K)^  which  again  leads  to  an  improved  error  estimate.  (Note  that 
Q*f^{K)  contains  extra  bubble  functions  compared  to  Qk{K),  but  condensing  out 
bubble  functions  is  relatively  cheap,  see  [24].) 

We  note  that  for  SQ,  MQl-2,4-6,  the  displacement  space  has  been  chosen  to  satisfy 

[Q't{k)f  C  v,(a:),  [qu,(K)]'  i.  Vi(A-). 
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This  ensures  that  all  the  displacement  spaces  have  the  same  0{h^)  asymptotic  ap- 
proximability  in  the  H^(f2)  norm  with  respect  to  the  h  version,  thus  facilitating  the 
easy  comparison  of  the  various  methods.  While  MQ3  does  not  have  this  property, 
it  is  useful  for  comparisons  with  SQ,  since  their  displacement  spaces  are  directly 
comparable. 

The  following  corollary  follows  from  Theorem  3.1  applied  to  parallelograms  and 
shows  that  the  above  choices  of  spaces  lead  to  mixed  methods  that  satisfy  a  common 
inf-sup  condition. 


Corollary  3.1  [20]  Let  the  mappings  Tk  be  affine,  i.e.  the  elements  are  all  parallel¬ 
ograms.  For  (Vat,  Wn)  defined  as  in  MQl-6  above, 


(3.6) 


inf  sup  r  - 

qeWN  veVjv  ||g||o  ||v||i 


(9,divv) 

n.ii  lu.ii"  =  pN>Ck  2, 


where  the  constant  C  is  independent  of  h  and  k. 


Corollary  3.1  shows  that  the  above  combinations  are  stable  with  respect  to  /i  as  /i  — >■  0. 
Moreover,  the  spaces  are  almost  stable  with  respect  to  A:  as  — >■  oo,  with  the  stability 
constant  behaving  like  Ck'~2,  (This  loss  is  unavoidable,  at  least  for  MQ5,  [4]). 


3.2  Numerical  Investigation  of  Inf-Sup  Condition 


For  a  given  mixed  method,  we  want  to  verify  computationally  that  inequality  (3.6) 
holds  and  since  it  is  a  lower  bound,  see  if  the  bound  is  approached.  We  proceed 
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as  in  [10].  Let  {(i>i}^i  be  a  basis  for  Wn{0,)  and  be  a  basis  for  Yn{0,). 

Then  q  G  ITiv(fi)  can  be  written  q  —  and  v  G  'Vn{Q,)  can  be  written 

V  =  E"=i  Then 

m  n 

{q,  div  v)  =  ajdiv  il^j)  =  div  ‘^■)cx  =  A'^Ba, 

i=i  j=i 

where  B  =  [(^,-,  div  V»j)]  is  m  x  n,  A  is  m  x  1  and  a  is  n  x  1.  Also 

m  m 

IkllS  =  (E  -'i*.  E  h)>'  =  •'’’DA, 

i—l  j—l 

where  D  =  ^i)]  is  m  x  m.  And 

n  n  n  n  n  n 

l|v||?  =  =  «^Ca, 

8=1  j—l  8=1  j=l  8  =  1  j  =  l 

where 


C  =  [{lpi,'>Pj)  +  +  (V’8-,2,'0j,2)] 


is  n  X  n.  The  discrete  inf —  sup  constant  in  (3.6)  can  then  be  characterized  as 


A^Ba 

mm  max  — - - j - r 

xeR-aeM"  (A^DA)2(arCa)2 


=  pv  >  0. 


Since  C  and  D  are  symmetric  positive  definite,  we  let 


2/  =  Dz A,  X  =  C^a. 


Then,  the  inf  —  sup  constant  p^  can  be  written 


y^T)-^BC-hx 

Pn  =  mm  max - - j-. 

3/eK”’a:eR”  (v^v)^(x^x)^ 
y:^0  ®#0  V»  W/  V  J 


28 


For  a  given  y,  the  maximum  over  x  occurs  when  y^D  2  BC' 


and  X  are  colinear,  or 


7C  2B^D  2y  =  X. 


Then 


and 


7y^D-lBC-iB^D-iy 

P;v  =  min  - i - i - ; - j- 

j/eR-  (y^y)2(72(yTi5-|Bc-iBTD-2y))l 


2  .  y^D-2BC-iB^D-ty 

p%  =  mm  ^ ^ 

yeR”" 
y^o 


y^y 


which  is  a  Rayleigh  quotient.  This  is  minimized  when 


D-iBC-‘B^D-5j,„  =  pi,y. 


or 

(3.7)  BC"^B^Mm  =  p%T>Um 

and  p^j^  is  the  smallest  positive  generalized  eigenvalue  of  (BC~^B^,D).  To  verify 
(3.6)  for  a  given  mixed  method  we  construct  the  matrices  B,  C  and  D  and  solve  (3.7) 
for  the  minimum  positive  eigenvalue. 

In  Figure  3.2,  we  show  the  results  of  numerically  estimating  the  inf-sup  constant  pN 
for  MQl  and  MQ5.  The  finite  element  spaces  are  defined  on  a  square  of  side  2,  with 
four  equal  square  elements  of  side  1  (with  natural  boundary  conditions).  The  0{k~^) 
deterioration  is  clearly  observed  for  both  elements,  showing  that  the  stability  of  MQl 
and  MQ5  is  similar.  For  the  given  mesh,  the  inf-sup  constant  for  MQ5  is  slightly  larger, 
about  1.09  times  the  size  of  the  inf-sup  constant  for  MQl.  Similar  computations  were 
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Figure  3.2:  The  inf-sup  constant  for  MQl  and  MQ5 

performed  in  [14]  where  they  considered  the  homogeneous  Dirichlet  problem.  The 
decay  rate  they  found  was  less  than  what  we  show  here.  When  comparing  the  actual 
magnitude,  their  inf-sup  constant  estimates  are  significantly  less  than  our  estimates, 
which  is  to  be  expected  since  the  displacement  space  for  the  Dirichlet  problem  is  a 
subspace  of  the  displacement  space  of  the  Neumann  problem. 

3.3  Asymptotic  Convergence  Rates 

is  the  minimal  conforming  space 
of  polynomials  of  degree  <  k.  This  means  that  the  asymptotic  convergence  rate  in 


We  note  that  in  SQ,  the  space  Yk{K)  = 
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terms  of  the  H^(0)  norm  for  SQ  is 

(3.8)  ||u  -  u„||,  <  C.  (u)  ||u||_^  _ 

assuming  the  solution  u  £  H’’  (0),  r  >  1.  This  is  the  asymptotic  rate  of  best  approx¬ 
imation  referred  to  in  Section  2.2.  A  method  that  produces  this  rate  for  all  u  G  [o,  ^ 
is  locking-free.  The  same  rate  will  also  be  observed  in  the  energy  norm  (2.5)  of  the 
error 

||u  -  ujvIIe,^  =  (2//||e  (u  -  Ujv)  ||o  +  A||div  (u  -  Uiv)||o) "  , 

which  is  equivalent  to  the  H^(Q)  norm  for  i/  bounded  away  from  1/2.  Also,  defining 
for  SQ 

(3.9)  pn  —  —A  div  uat, 
we  see  that 

(3.10)  Up  -  p«||„  <  c,  (p)  ||u||^ . 

The  upper  bound  in  (3. 8), (3. 10)  represents  the  optimal  asymptotic  convergence 
rate  for  polynomials  of  degree  k.  Unfortunately,  both  constants  Ci  and  C2  blow  up  as 
1/  — )■  I  (A  — 00),  so  that  this  rate  is  not  uniform  if  u  is  allowed  to  vary  in  [0,  |).  The 
following  theorem  gives  the  best  possible  estimate  independent  of  u  for  the  special 
case  of  a  uniform  rectangular  mesh. 

Theorem  3.2  [1]  Let  u  ^  H’’  (O),  r  >  1,  be  the  solution  of  (2.3)  and  let  un  be  the 
solution  obtained  by  SQ,  using  parallelogram  elements.  Then  there  exists  a  constant 
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C  independent  ofu,  h,  k  and  u  €  [0,^)  such  that 

(3.11)  ||u  -  ujvlli  <  ^ 

where  a  =  1  for  fc  =  1, 2  and  a  —  2  for  k  >  2.  Moreover,  this  is  the  sharpest  estimate 
that  is  uniform  in  v  for  v  G  [0,  ^). 

The  same  estimate  in  equation  (3.11)  also  holds  in  the  energy  norm  defined  by  (2.5). 
For  results  on  other  meshes,  see  [1]. 

We  see  from  the  above  theorem  that  if  the  h  version  is  used,  the  rate  deteriorates 
by  O  {hf)  as  u  ~(0  (h)  if  A;  =  1,2).  In  terms  of  Definition  2.2  h-refinement  for  SQ 
has  0{h~'^)  locking  for  A:  >  2  and  0  {h~^)  locking  for  A:  =  1,2  in  the  displacement. 
For  the  p  version,  there  is  no  deterioration  in  the  asymptotic  convergence  rate,  i.e., 
p-refinement  for  SQ  is  locking-free  in  the  displacements  in  the  sense  of  Definition  2.2. 

Remark  3.2  The  statement  that  there  is  no  locking  with  the  p  version  for  displace¬ 
ments  is  somewhat  misleading.  Even  though  there  is  no  deterioration  in  the  asymp¬ 
totic  rate  of  convergence,  the  actual  observed  error  does  deteriorate,  particularly  with 
curved  elements.  In  other  words,  the  locking  ratios  defined  in  Definition  2.1  increase. 
See  [22]  and  Chapter  5. 

Let  us  mention  that  if  in  SQ  we  take  Vk{K)  =  |^(5fc(.^)j  instead,  then  a  =  1  for 
all  k  in  (3.11),  so  that  this  choice  has  O  {h)  locking  for  h-refinement  for  all  k. 

Turning  now  to  the  mixed  methods,  we  have  the  following  theorems,  which  follow 


from  Theorem  3.1. 
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Theorem  3.3  [20]  Let  (u,p)  the  solution  of  (2.7)-(2.8)  be  in  H'’ (fi)  x  H’’~^  (0), 
r  >  1.  Let  (ujVjPiv)  be  the  solution  of  MQ5  solving  equations  (2.9)-(2.10)  with  paral¬ 
lelogram  elements.  Then  for  any  e  >  0,  there  exists  a  constant  independent  of  u, 
h,  k  and  v  G  [0,  |)  such  that 

(3.12)  ||u  -  u„||,  <  (||u||^  +  IIpILJ  , 

(3.13)  Up  -  p„||„  <  (||u||^  +  iipII,^,)  . 

Theorem  3.4  [20]  Let  (u,p)  be  as  in  Theorem  3.3  and  (uAr,pjv)  be  the  solution  of 
MQl-4,6  solving  equations  (2.9)-(2.10)  with  parallelogram  elements.  Then  for  any 
€  >  0,  there  exists  a  constant  independent  of  u,  h,  k  and  u  G  [0,  |)  such  that 

(3.14)  ||u  -  u«||,  <  (Hull,  +  IIpILJ  , 

(3.15)  Up  -  PA,||„  <  . 

From  the  above,  we  observe  that  in  terms  of  the  p  version,  all  the  mixed  methods 
are  (like  SQ)  asymptotically  free  of  locking  in  displacement  (except  for  an  O  (k^) 
factor,  which  comes  in  due  to  a  technicality  in  the  proof).  In  terms  of  the  h  version, 
MQ5  shows  an  O  {h~^)  loss  in  the  displacement  rate  of  convergence  when  compared  to 
the  optimal  estimate  (3.8)  (for  r  >  A;  +  1).  Similarly,  MQ3  shows  an  0{h~^)  loss  in  the 
displacement  rate,  since  with  elements,  we  would  expect  0(/j*+^)  convergence  in 
the  displacements  for  smooth  solutions.  However,  these  losses  are  not  manifestations 
of  locking,  since  the  loss  is  observed  for  any  u.  Rather,  this  is  an  approximability 
problem,  since  the  space  for  p^  should  be  of  one  degree  lower  than  that  for  Ujv  (and 
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not  two  degrees  lower).  Hence,  all  these  methods  are  free  of  locking  in  displacement 
for  h-refinement,  and,  in  addition,  MQ1,2,4,6  are  correctly  balanced. 

Turning  to  the  stresses,  for  the  h  version,  the  conclusions  for  locking  are  the  same 
as  those  for  displacements.  For  the  p  version,  if  we  compare  (3.13),  (3.15)  with  the 
optimal  rate  in  (3.10),  we  see  a  loss  of  Hence,  we  conclude  that  p-refinement 

gives  rise  to  stress  locking  of  order  0(A:~^2+^))  for  all  the  above  mixed  methods.  As 
we  see  in  the  next  section,  this  is  a  rather  mild  amount  of  locking,  and  does  not  show 
up  in  the  model  problems  for  the  practical  ranges  of  k  used  by  us,  giving  a  method 
that  is  essentially  locking-free. 

Let  us  mention  that  (3. 12), (3. 14)  will  not  hold  with  ||u  — UAr||i  replaced  by  the  error 
||u  — UatHe',;/  defined  by  (2.5).  The  correct  analog  of  the  “energy  norm  error”  for  mixed 
methods  is  given  in  (2.11).  With  this  definition  of  the  energy  norm,  (3. 12), (3. 14)  will 
hold  once  again. 

Remark  3.3  With  proper  mesh-degree  selection^  exponential  rates  of  convergence 
may  he  obtained  using  the  hp  version.  The  estimates  in  Theorems  3.2  -  3.4  would 
then  be  improved,  reflecting  these  exponential  rates.  Exponential  convergence  requires 
highly  refined  meshes  in  parts  of  the  domain  (e.g.  at  corners).  Such  meshes  may  he 
constructed  even  if  the  elements  are  all  parallelograms  by  using  “hanging  nodes.  ”  See 
[11]  for  a  discussion  of  a  fully  adaptive  hp  method  (with  hanging  nodes)  for  the  Stokes 
problem,  based  on  MQ].  Also,  see  [17]  where  the  hp-FEM  is  constructed  by  using  both 
triangular  and  parallelogram  elements  to  give  such  exponential  rates  of  convergence 
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on  polygonal  domains. 

Remark  3.4  The  mixed  methods  described  here  may  be  generalized  to  three-dimensional 
parallelepiped  elements,  see  [20].  The  best  stability  constant  in  the  3-d  case  is  Ck~^ . 

Remark  3.5  In  [20],  Theorems  3.1,  3.3  and  3.]  above  have  actually  been  established 
for  the  more  restrictive  case  of  homogeneous  Dirichlet  boundary  conditions.  In  that 
case,  the  spaces  for  the  displacements  are  further  restricted,  so  we  can  expect  the  lock¬ 
ing  to  be  more  severe.  In  subsequent  sections,  it  is  shown  that  the  MFEMs  perform 
much  better  than  STEM  with  respect  to  locking,  using  test  problems  involving  Neu¬ 
mann  (traction)  boundary  conditions.  For  the  case  of  a  clamped  boundary,  we  expect 
the  difference  in  performance  to  be  even  more  significant. 

3.4  Numerical  Experiments 

In  this  section,  we  investigate  the  computational  performance  of  the  methods  SQ, 
MQl-5,  both  with  h-  and  p-refinement.  Our  results  support  the  theoretical  evidence 
in  Section  3.3,  i.e.  that  only  appropriately  balanced  mixed  methods  provide  near- 
optimal  performance  both  for  displacements  and  stresses  calculated  either  by  h-  or 
p-refinement.  Although  we  did  not  perform  computations  with  MQ6,  we  expect  its 
numerical  behavior  to  be  similar  to  MQl-4,  as  suggested  by  Theorems  3.1  and  3.4. 
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3.4.1  Model  Problem 


In  our  numerical  experiments,  we  use  a  model  problem  that  we  refer  to  as  SING, 
having  reference  to  the  movable  singularity  in  the  displacement  u.  For  this  problem, 
we  solve  (2.4)  or  (2.9)  -  (2.10)  on  0  =  {(a^,?/)  :  — 1  <  a;  <  1,  — 1  <  y  <  1},  with  f  =  0 
and  tractions  g  specified  on  the  boundary  by 

8(a;-a;o)(y-yo)[(y~yo)^-(a;-a;o)^ 

^ 

4  (x  -  XqY  [3  (y  -  yof  -  {x  -  Xq)^] 

(x-xo)^  +  (y-yo)^]^ 

8  (x  -  xo)  (y  -  yo)  [(y  -  yof  -  {x  -  xo)^] 

g2{x,y)  =  - - r - + 

(x  -  Xo)  +  (y  -  yo) 


4  (y  -  yof  \{y  ~  yof  -  3  (x  -  xo)^ 

- F - - - t3 - ^”2- 

[(x  -  xo)^  +  (y  -  yof\ 

Here  (ni,n2)  is  the  outward  unit  normal  on  dfl.  The  exact  solution  is  given  by  [15] 


(3.16) 


(3.17) 


ui  (x,y) 


«2(a;,y)  = 


(x  -  Xo)  (A  +  2y)  (x  -  xo)^  -  A  (y  -  yo)^ 
y  (A  +  y)  [(x  -  xo)^  +  (y  -  yo)^]^ 

(y  -  yo)  [a  (x  -  Xo)^  -  (A  +  2y)  (y  -  yof 
y  (A  +  y)  [(x  -  xo)^  +  (y  -  yo)^]^ 


(from  which  the  stresses  can  be  calculated).  In  the  above,  A  and  y  are  defined  in 
terms  of  E  and  i/  by  (2.2).  We  take  =  1  and  let  v  vary. 

As  seen  in  (3.16)  -  (3.17)  the  solution  has  a  movable  singularity  at  (xo,  yo)-  By 
taking  (xo,yo)  =  (-2,-2),  we  get  the  case  of  a  smooth  solution.  We  also  consider 
the  case  (xo,  yo)  =  (“l-l,  —l-l)?  to  see  how  well  these  methods  handle  a  near-singular 


solution. 
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3.4.2  Computational  Considerations 

Implementation,  of  our  mixed  method  is  relatively  easy.  Since  no  inter-element  con¬ 
tinuity  is  imposed  on  pjv,  we  simply  solve  (2.10)  locally  over  each  element  for  pjv  in 
terms  of  u^v,  and  then  substitute  this  in  (2.9).  This  gives  a  positive  definite  system 
for  Uiv,  which  can  then  be  solved  by  a  direct  solver.  For  the  range  of  i/’s  tested 
(0.3  <  u  <  0.4999),  round-off  error  was  not  a  problem. 

The  MFEMs  do  cost  more  than  the  SFEMs  to  compute.  This  increase  is  due  to 
the  pressure  variable,  which  increases  the  size  of  the  element  stiffness  matrix.  Eor 
example,  for  MQl  the  increased  cost  ranges  from  36%  when  k  =  2  to  75%  when 
k  =  10,  up  to  125%  as  fc  — >■  oo.  However,  the  overall  cost  of  the  two  methods  is  more 
difficult  to  compare.  Since  the  pressures  are  discontinuous,  they  can  be  condensed 
out  on  each  element.  Similarly,  if  a  frontal  solver  is  used[24],  the  pressure  degrees 
of  freedom  do  not  contribute  to  the  size  of  the  front.  Another  factor  is  the  size  of 
the  problem.  For  problems  with  few  elements,  the  MFEM  cost  may  be  significant. 
But  when  solving  large  problems,  the  cost  of  solving  the  global  system,  which  is  the 
same  for  MFEM  and  SEEM,  may  dominate.  On  a  parallel  machine,  cost  comparisons 
become  even  less  obvious. 

Eor  these  reasons,  rather  than  comparing  the  error  obtained  via  different  methods 
with  N  the  sum  of  dim(V^)  and  dim(lFjv),  we  compare  it  with  n  =  dim(VAr),  which 
represents  the  total  number  of  variables  in  the  global  system  to  be  solved.  Of  course, 
our  choice  n  is  just  one  possible  (incomplete)  measure  of  the  actual  cost,  which  is  very 
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implementation  dependent.  The  SFEM  will  definitely  be  cheaper  when  locking  is  not 
an  issue.  The  MFEMs  are  more  robust,  but  are  cheaper  only  when  1/  is  sufficiently 
close  to  0.5. 

The  mesh  Ch  is  a  uniform  partition  of  0  into  square  elements  of  side  h.  Then  it  is 
easily  seen  that  the  number  of  degrees  of  freedom  N  satisfies 

N{h,k)  ^  Ch-‘^k^, 

so  that  the  asymptotic  rates  of  Section  3.2  could  be  re-written  (for  comparison)  in 
terms  of  iV,  using  the  equivalence  h  for  h-refinement  {k  fixed)  and  k~^  ~ 

for  p- refinement  {h  fixed). 

3.4.3  Computational  Convergence  Results 

Comparison  of  SQ  and  MQl  -  h  version 

Figures  3.3  -  3.5  show  the  results  of  h-refinement  with  SQ  and  MQl  (smooth  solution). 
The  number  of  elements  ranges  from  1  to  64.  We  have  plotted  the  percentage  relative 
errors  in  the  energy  norm  of  the  displacement  (defined  by  (2.5)  and  (2.11)),  and  the 
norm  of  the  sum  of  the  normal  stresses  (SNS),  vs  n,  the  number  of  degrees  of  freedom. 
(Recall  that  the  energy  norm  satisfies  the  same  error  estimate  as  the  norm  of  the 
displacement.)  On  these  h-refinement  plots  “rate”  is  the  exponent  of  h,  i.e.  the  rate 
of  convergence. 

For  SQ,  the  energy  norm  error  shows  good  correlation  with  the  asymptotic  rates 


predicted  in  Equation  3.8 


the  optimal  rate  of  O  for  =  2, 3, 4  is  observed 
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Figure  3.3;  SQ  vs  MQl,  /j-refinement,  Smooth  solution,  k  =  2. 


Figure  3.4;  SQ  vs  MQl,  /j-refinement.  Smooth  solution,  k  =  3. 

when  V  —  0.3,  while  the  “locked”  rate  of  O  for  k  =  2, 3  is  observed  when 

V  =  0.4999.  When  A;  =  4,  we  observe  the  predicted  locking  of  order  0{h~'^)  on  this 

model  problem.  The  energy  norm  locking  ratios  ^  {Le)  in  Table  3.4.3  emphasize  that 

the  locking  ratio  increases  gradually  as  N  increases  {h  decreases).  In  contrast,  MQl 
^The  locking  ratio  if; (0.4999,  iV)  is  the  observed  error  a,t  i/  —  0.4999  divided  by  the  observed 


error  ai  v  =  0.3. 
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exhibits  the  asymptotic  rate  of  best  approximation  in  energy  norm  for  k  =  2,3  and  4 
,  and  is  locking-free.  The  locking  ratios  (not  shown)  remain  bounded  by  1.15. 


Figure  3.5:  SQ  vs  MQl,  h-refinement,  Smooth  solution,  k  =  4. 


For  SQ,  the  deterioration  in  SNS  is  worse,  since  there  is  no  observed  convergence 
for  fc  =  2  when  u  =  0.4999.  The  error  is  just  beginning  to  decrease  for  A:  =  3  in 
the  range  of  h  we  use.  When  k  =  4,  SQ  begins  to  converge  for  stresses  with  a  rate 
of  O  The  stress  locking  ratios  ^  (Ls)  in  Table  3.4.3  show  dramatic  increase 

with  N.  On  the  other  hand,  we  see  that  MQl  produces  the  same  near-optimal  results 
for  both  u  =  0.3  and  0.4999.  The  stress  locking  ratios  (not  shown  for  MQl)  remain 
bounded  by  1.04. 


N 


N 

13 

39 

m 

127 

189 

263 

349 

447 

557 

679 

Lb  (0.4999,  iV) 

1.0 

1.7 

4.0 

6.2 

7.0 

£s(0.4999,iV) 

1.2 

10.3 

20.0 

33.6 

50.9 

71.7 

95.8 

122.9 

152.5 

184.5 

k  =  3 

N 

21 

63 

125 

207 

309 

431 

573 

735 

917 

1119 

(0.4999,  AT) 

1.5 

1.8 

2.6 

3.8 

5.4 

8.9 

10.5 

11.8 

13.0 

is  (0.4999,  iV) 

11.9 

23.4 

52.1 

106.4 

188.4 

290.2 

396.7 

494.8 

578.7 

11 

N 

m 

95 

191 

319 

479 

671 

895 

1151 

1439 

1759 

iB(0.4999,Ar) 

1.7 

3.6 

4.9 

6.2 

7.3 

9.3 

10.2 

is  (0.4999,  AT) 

9.3 

12.3 

40.4 

100.0 

189.6 

293.7 

394.1 

481.2 

554.7 

618.0 

Table  3.1:  Locking  Ratios,  SQ,  /i-refinement,  k  =  2,  3,  4. 
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Comparison  of  SQ  and  MQl  -  p  version 

Next,  in  Figures  3.6  and  3.7  we  see  the  perfornaance  of  SQ  and  MQl  with  p  -refinement 
on  a  four  element  mesh.  We  observe  no  locking  with  either  method  where  the  displace¬ 
ment  is  concerned,  since  the  curves  for  u  =  0.3  and  u  =  0.4999  are  very  close.  (Here, 
and  for  all  mixed  method  calculations,  the  “energy”  norm  used  is  the  one  defined  in 
(2.11).)  If  one  looks  at  the  stresses  for  the  smooth  solution,  there  is  a  marked  de¬ 
crease  in  accuracy  with  SQ  when  we  go  from  u  =  0.3  to  =  0.4999.  This  decrease  in 
accuracy  is  even  more  apparent  for  the  unsmooth  solution,  see  Table  3.4.3 


2  3 


Degrees  of  Freedom  Degrees  of  Freedom 

Figure  3.7:  SQvsMQl,  p-refinement.  Unsmooth  solution. 

On  the  other  hand,  MQl  shows  no  change  as  i/  — >  0.5,  either  in  the  smooth  or 
non-smooth  case.  Hence,  no  locking  can  be  seen  in  these  results  for  MQl,  not  even 

the  O  predicted  by  Theorem  3.4.  The  results  with  the  other  mixed  methods 

®The  locking  ratio  2/5(0.4999,  iV)  is  the  observed  error  at  t/  =  0.4999  divided  by  the  observed 


error  at  i/  =  0.3. 
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Table  3.2:  Locking  Ratios,  SQ,  p-refinement,  Smooth,  Unsmooth. 


are  similar,  showing  that  stresses  are  accurately  recovered  by  the  p  version  mixed 
method.  (Let  us  remark,  however,  that  even  though  the  O  loss  is  not  seen 

for  this  problem,  it  could  appear  in  other  problems,  depending  upon  the  nature  of  the 
exact  solution.) 

Comparison  of  MQl-5  -  h-  and  p-refinement 

Next,  we  compare  MQl-5,  by  seeing  how  well  they  perform  with  h-refinement.  As 
mentioned  in  the  introduction,  for  an  hp  code,  optimality  with  respect  to  h  as  well  is 
highly  desirable.  In  Figures  3.8  and  3.9,  we  plot  the  displacement  and  stress  errors 
of  MQl-5  when  the  solution  is  smooth,  for  k  =  2  and  4.  These  have  been  plotted  for 
u  =  0.4999  (the  figures  do  not  change  appreciably  for  i/  closer  to  |  or  for  u  =  0.3). 
We  see  clearly  the  superiority  of  MQl-4,  which  have  0{h^)  convergence  in  both 


Figure  3.8:  MQl-5,  /j-refinement,  Smooth  solution,  k  =  2. 


Figure  3.9:  MQl-5,  /i-refinement,  Smooth  solution,  k  =  4. 


energy  norm  and  SNS,  which  is  the  optimal  rate  for  MQ1,2,4,5.  The  reason  MQ5 
gives  a  sub-optimal  rate  (0(/i^“^))  is  not  due  to  locking,  since  the  rate  is  the  same  for 
u  =  0.3  and  i/  =  0.4999.  Rather,  it  is  due  to  the  polynomial  degree  nsed  for  pj\j  being 
2  degrees  lower  than  that  for  Ujv,  i.e.  it  is  an  approximability  issne.  One  could,  of 


course,  compare  MQl-4  with  degree  k  to  MQ5  with  degree  1,  but  then  the  optimal 
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rate  expected  for  the  displacements  from  MQ5  would  be  one  order  higher,  and  would 
still  not  be  attained.  Note  that  the  optimal  rate  for  MQ3  is  convergence.  It 

fails  to  achieve  this  rate  for  the  same  reason  MQ5  does,  because  the  pressure  spaces 
is  3  degrees  lower  than  the  displacement  space.  This  emphasizes  the  importance 
of  considering  both  stability  and  approximability  when  constructing  mixed  method 


spaces. 


Figure  3.10;  MQT5,  p-refinement,  Smooth  solution. 


For  completeness,  we  also  compare  MQl-5  under  p-refinement  (for  four  elements). 


These  results  (Figures  3.10  and  3.11)  show  the  mixed  methods  perform  about  equally 


well.  Since  we  have  not  taken  dim(VF7v)  into  accormt  in  the  comparisons,  the  method 


with  the  minimal  degrees  of  freedom  for  i.e.  MQl,  may  be  preferable  overall. 


%  Rel  Error  -  Energy  Norm 


Chapter  4 


Straight-Sided  Triangular  Elements 


From  the  results  in  [26],  it  follows  that  the  inf  —  sup  constant  for  the  p  version  SFEM 
on  triangular  elements  is  0{k~'*)  for  some  unknown  7.  Many  commercial  codes  rely 
heavily  on  automatic  mesh  generation  using  triangular  elements.  For  these  codes  to 
have  reliable  hp  capability,  triangular  elements  with  known  stability  characteristics 
are  necessary.  In  this  chapter,  we  look  at  mixed  method  triangular  elements  for  which 
a  better  characterization  of  stability  can  be  made.  These  mixed  triangular  elements 
are  based  on  conditions  (A1)-(A6)  outlined  in  Chapter  3.  The  pair  of  spaces  that 
was  analyzed  by  Schwab  and  Suri  [17]  satisfy  these  conditions.  They  showed  that  the 
inf  — sup  constant  is  no  worse  than  Ck~^.  We  provide  evidence  that  the  inf  — sup 
constant  for  this  mixed  triangular  element  is  better  by  computing  it  for  a  single 
element  for  a  practical  range  of  k,  and  then  invoking  Theorem  3.1.  We  also  propose 
and  investigate  an  optimized  version  of  this  element.  We  then  perform  experiments  on 
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triangular  meshes  with  both  SFEM  and  our  mixed  triangular  elements  and  compare 
their  performance  for  a  model  problem. 

4.1  Spaces  and  Stability 

As  in  section  3.1,  we  use  the  same  definition  of  Vn  and  Wn,  with  the  exception  that 
now  the  reference  element  K  =  {(x,  y):0<x<l,0<y<l  —  a;}isa  triangle.  We 
use  conditions  (A1)-(A6)  (see  section  3.1)  to  select  displacement  and  pressure  spaces 
for  two  mixed  triangular  elements. 

We  identify  the  first  triangular  element  as  MTl,  and  begin  by  letting  Y k{k)  = 
Vk{k)  .  Then  (Al)  is  satisfied  and  (A2)  requires  Wk{k)  C  Vk-2{k),  so  we  let 
Wk{K)  =  Vk-2{K).  Condition  (A3)  is  verified  in  [17]  for  these  spaces  with  7  =  3. 
Condition  (A4)  is  satisfied  if  V2{k)  C  Vk{k).  Finally,  (A5)  is  satisfied  and  (A6) 
is  satisfied  with  n  =  1.  We  point  out  that  n  >  1  in  (A6)  implies  less  than  optimal 
approximation.  Thus  we  define 

MTl:  Vk{k)  =  [Vk{k)]\  Wk{k)  =  Vk-2{k), 

and  expect  MTl  to  have  inf-sup  constant  no  worse  than  0{k~^)  and  less  than  optimal 
approximability. 

Another  approach  is  to  set  Wk{K)  =  Vk-i{K)  and  determine  the  smallest  space 
Vfc(jA)  that  satisfies  conditions  (Al)  -  (A6).  The  lowest  order  triangular  bubble  func¬ 
tion  e  Vz{k),  and  '^q  €  Vk-i{k)  =  Wk{k)  we  have  ^  'Pk-2{k).  Then  (A2) 
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requires  that  {bubble  functions  in  C  Vk{K).  This  leads  to  the  following 

mixed  triangular  element 

MT2:  Vfc(A0  =  [VkiK)  U  {Vk+iiK)  bubble  functions}]^ 

Wkik)  =  VkMK). 

The  analysis  of  local  stability  for  MTl  also  applies  to  MT2,  so  the  inf-sup  constant 
will  be  no  worse  than  0{k~^).  Note  that  in  contrast  to  MTl,  the  polynomial  degree 
of  the  displacement  and  pressure  spaces  for  MT2  differs  by  only  1,  so  MT2  should 
have  optimal  approximability. 

For  SFEM,  we  define 

ST:  Vk{k)  =  Vk(k). 

This  is  equivalent  to  SQ,  since  for  a  given  k,  both  SQ  and  ST  are  the  minimal  space  of 
polynomials  of  degree  k  that  produces  C°  continuity  on  a  regular  mesh.  As  mentioned 
earlier,  the  inf-sup  constant  for  ST  is  0{k~'^)  for  some  unknown  7  [26]. 

To  characterize  the  stability  of  MTl  and  MT2,  we  would  like  to  use  Theorem  3.1. 
Both  MTl  and  MT2  satisfy  (Al)  and  (A2).  They  also  satisfy  (A3)  with  7  =  3 
(computations  indicate  this  is  pessimistic).  Condition  (A4)  is  met  by  both  MTl  and 
MT2  (see  [6]).  To  gain  a  practical  description  of  the  inf-sup  constant  for  these  mixed 
methods  we  turn  to  the  computational  approach  of  the  next  section. 
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4.2  Numerical  Investigation  of  Inf-Sup  Condition 

Recall  from  Chapter  3  that  the  global  inf-sup  constant  in  Theorem  3.1  can  be  es¬ 
tablished  if  either  (Al)-(A4)  hold  or  if  (LS),(A4)  hold.  Just  like  we  estimated  the 
global  inf-sup  constant  by  an  eigenvalue  computation  in  section  3.2,  so  also  we  can 
estimate  the  local  inf-sup  constant  in  (LS)  instead.  Once  7  in  (LS)  has  been  estimated 
computationally,  we  can  expect  the  estimate  in  Theorem  3.1  to  hold  for  the  range  of 
k  tested,  since  (A4)  holds. 

We  therefore  estimate,  for  both  MTl  and  MT2,  the  value  of  7  in  (LS)  over  the 
reference  triangle  K.  The  results  are  given  in  Figure  4.2.  We  see  that  for  A:  <  11  the 
slope  is  —0.3  and  for  A;  >  11  the  slope  is  —1.1.  This  significant  shift  in  the  slope  may 
indicate  that  the  A:  >  11  region  reflects  the  asymptotic  rate.  This  suggests  that  the 
inf-sup  constant  for  MTl  is  about  0{k~^). 

Most  commercial  codes  limit  p-refinement,  using  an  upper  bound  on  polynomial 
degree  that  is  usually  much  less  than  20.  For  practical  purposes  therefore,  the  0{k~^) 
computational  estimation  we  have  obtained  will  govern  the  globally  observed  inf-sup 
constant. 

Remark  4.1  Let  us  remark  that  ST  could  also  be  put  in  the  form  of  a  mixed  method, 
i.e.  in  the  form  (2.7)-(2.8),  by  taking  Wn  =  divfV^),  which  gives  [19]  Wk{K)  = 
Vk-i{K).  We  could  then  attempt  to  numerically  estimate  the  global  inf-sup  constant 
by  the  local  one  as  done  here,  since  Theorem  3.1  will  again  hold.  Unfortunately,  the 
local  inf-sup  constant  is  0  in  this  case,  even  though  the  global  one  is  positive.  This 
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Figure  4.1:  The  local  inf-sup  condition  for  MTl  and  MT2. 
discrepancy  occurs  since  condition  (A2)  is  not  satisfied  here. 


4.3  Asymptotic  Convergence  Rates 

Recall  that  the  space  Vfc(^)  =  used  in  ST  is  the  minimal  conforming  space 

of  polynomials  of  degree  <  k.  As  in  Section  3.2,  this  means  that  in  the  absence  of 
locking  {u  fixed  and  not  close  to  1/2),  the  asymptotic  rate  of  best  approximation  in 
terms  of  the  norm  for  ST  is  given  by  (3.8),  assuming  the  solution  u  G  H''  (0), 

r  >  1.  The  same  rate  will  also  be  observed  in  the  energy  norm  (2.5).  Also,  the 
inequality  given  in  (3.10)  holds  for  the  pressure.  The  upper  bormds  in  (3. 8), (3. 10) 
represent  the  optimal  asymptotic  convergence  rate  for  polynomials  of  degree  k. 

Keeping  h  fixed  in  (3.8),  we  get  the  rate  for  the  displacement  error  for  ST  un- 
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der  p-refinement  for  any  v  G  [0,  |).  Hence  p-refinement  is  free  of  locking  for  the 
displacements  (but  not  necessarily  the  pressures). 

Under  /i-refinement  with  fc  <  3,  0(h)  energy  norm  locking  is  the  best  rate  in¬ 
dependent  of  f.  However,  for  A:  >  4,  the  rate  of  convergence  in  displacement  and 
pressure  for  ST  is  independent  of  z/,  where  the  solution  u  €  H^. 

For  the  mixed  methods,  we  have  the  following  result  for  the  case  of  straight  sided 
triangles. 

Theorem  4.1  Let  V7v(H),  VUv(O)  (k  >  2)  satisfy  (LS),(A4).  Then  for  any  e  >  0, 
there  exists  a  constant  independent  ofu,h,k  and  u  G  [0,  |)  such  that 

||u  -  Ua^IIi  <  +  Iblls-i), 

lb-Piv||o  <  +  Ibll^-i), 

where 


I  =  min{s  —  1,A:  —  n},  n  =  0  for  MT2,n  =  1  for  MTl. 
and  7  is  the  exponent  in  (3.5). 

Proof:  We  note  that  since  (A4)  is  satisfied.  Theorem  3.1  holds.  The  proof  then  follows 
the  same  line  as  that  of  Theorem  5.2  in  [20]  for  parallelograms.  □ 

As  mentioned  previously,  we  have  characterized  the  constant  7  computationally, 
which  will  be  the  same  over  all  the  (affine-equivalent)  elements. 


4.4  Numerical  Results 
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In  this  section  we  examine  the  performance  of  ST  and  MTl  on  model  problem  SING 
(see  Section  3.4.1).  We  are  interested  in  both  h-  and  p-refinement  and  also  consider 
both  smooth  and  near-singular  solutions.  We  compare  the  results  to  the  error  bounds 
given  above. 

4.4.1  /t- Refinement 

We  begin  with  ST  using  h-refinement  on  both  a  smooth  and  a  near-singular  solution. 
We  see  in  Figure  4.2  that  when  the  solution  is  smooth  and  k  =  2,  energy  norm  locking 
of  0(h)  is  observed  as  z/  — This  order  of  locking  is  predicted  by  (3.8)  for  A;  =  3  as 
well,  although  it  does  not  appear  for  this  model  problem.  Consistent  with  theory,  we 
see  no  energy  norm  locking  when  A:  =  4. 


Degrees  of  Freedom  —  h-refinement  Degrees  of  Freedom  —  h-refinement 

Figure  4.2;  ST,  h-refinement.  Smooth  Solution. 


The  SNS  shows  no  convergence  for  A;  =  2  (i.e.  O(h^)  locking),  and  0(h)  locking 
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Figure  4.3:  ST,  /i-refinement,  Unsmooth  Solution. 


for  k  =  3.  When  fc  =  4,  the  SNS  rate  of  convergence  is  unaffected  as  but  we 

see  a  locking  ratio  of  about  4  that  increases  gradually  as  h  decreases,  see  Table  4.4.1. 

When  the  solution  is  near  singular,  we  see  in  Figure  4.3  the  rate  of  convergence  in 
energy  norm  does  not  change  as  u  This  is  consistent  with  (3.8),  since  the  rate 

depends  both  on  k  and  the  regularity  of  the  solution.  For  near  singular  solutions,  the 
effective  r  is  quite  small,  so  the  rate  of  convergence  fj,  =  min{A;,r  —  1}  is  dominated 
by  r.  The  behavior  for  the  pressure  is  qualitatively  the  same  for  the  near-singular 
solution  as  for  the  smooth. 

Next  we  use  MTl  under  /i-refinement  to  solve  SING  with  a  smooth  solution  first, 
and  then  a  near-singular  one.  In  Figures  4.4  and  4.5  we  see  that  the  performance  of 
MTl  is  unaffected  by  the  value  of  i/.  Even  the  pressure  approximation  is  unaffected 
as  z/  As  with  ST,  the  regularity  of  the  near-singular  solution  dominates  the  rate 

of  convergence  in  Figure  4.5,  so  even  for  k  =  4,  the  observed  energy  norm  rate  of 


k  =  2 


5  47 


159  239  335 


is  (0.4999,  iV) 

1.37 

2.97 

5.39 

8.17 

11.35 

14.94 

18.95 

CO 

II 

N 

29 

95 

197 

335 

509 

719 

965 

1247 

Ts  (0.4999,  iV) 

1.24 

2.71 

4.40 

5.56 

6.74 

8.11 

9.75 

11.66 

k  =  4 

N 

47 

159 

335 

575 

879 

1247 

1679 

2175 

Ls{0A999,N) 

1.34 

2.36 

3.21 

3.61 

3.86 

4.03 

Table  4.1:  SNS  Locking  Ratios,  ST,  n-rehnement,  Smooth  Solution,  k 
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convergence  never  exceeds  0{h^). 

Note  that  similar  to  MQ5,  MTl  under  /i-refinement  is  locking  free,  yet  the  rate 
of  convergence  in  both  energy  norm  and  SNS  is  0{h)  worse  than  the  asymptotic  rate 
of  best  approximation.  This  is  due  to  the  polynomial  degree  for  pjv  being  2  degrees 
lower  than  that  for  u^r  (rather  than  1  as  in  MT2). 


Degrees  of  Freedom  Degrees  of  Freedom 


Figure  4.5:  MTl,  h-refinement.  Unsmooth  Solution. 

4.4.2  ^-Refinement 

We  see  in  Figures  4.6  and  4.7,  that  under  p-refinement,  ST  performs  very  well  on 
our  model  problem.  The  displacement  approximation  is  essentially  the  same  for  both 
1/  =  0.3  and  u  =  0.4999.  Only  the  pressures  show  a  slight  difference,  with  a  nearly 
constant  locking  ratio  of  about  1.4  for  the  near-singular  solution  on  8  elements. 

Even  this  small  shift  is  avoided  by  MTl.  We  see  in  Figures  4.8  and  4.9  that  the 
value  of  u  does  not  affect  the  performance  of  MTl. 
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Figure  4.6:  ST,  p-refinement,  Smooth  Solution. 


Degrees  of  Freedom  —  k  =  2  ...  10 


Degrees  of  Freedom  —  k  =  2  ...  10 


Figure  4.7:  ST,  p-refinement.  Unsmooth  Solution. 


In  summary,  the  numerical  performance  of  ST  and  MTl  on  our  model  problem 
is  consistent  with  available  theory.  ST  locks  in  the  energy  norm  for  low  k  under  h- 
refinement,  but  is  locking  free  if  A:  >  4.  As  predicted,  for  the  p  version,  it  is  locking-free 
in  the  energy  norm  but  (somewhat  surprisingly)  almost  locking-free  in  the  pressures 
as  well.  Pressure  locking  is  more  pronounced  in  the  h  version.  Here,  no  convergence 
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Figure  4.8;  MTl,  p-refinement,  Smooth  Solution. 


Figure  4.9:  MTl,  p-refinement,  Unsmooth  Solution. 


is  observed  for  k  =  2  and  0{h)  locking  for  fc  =  3,  for  both  the  smooth  and  near¬ 
singular  solutions.  When  A:  =  4  no  pressure  locking  occurs,  although  the  locking  ratio 
may  be  significant  for  smooth  problems.  The  experiments  indicate  that  MTl  under 
h-refinement  is  locking  free  for  any  reasonable  k. 

In  light  of  these  experiments,  the  lower  bound  on  the  inf —  sup  constant  for  MTl 
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obtained  analytically  in  [17]  seems  pessimistic.  The  inf-sup  estimate  of  k~^'^  for  A;  <  11 
suggested  by  the  computations  of  Section  4.2  is  reasonable  since  MTl  shows  no  SNS 
locking  under  p-refinement  in  that  range.  To  test  the  apparent  shift  suggested  by 
Figure  4.2,  the  p-refinement  experiments  of  Figures  4.8  and  4.9  on  MTl  would  need 


to  be  run  with  higher  values  of  k. 


Chapter  5 


Curvilinear  Elements 


Since  curved  elements  are  essential  for  p  version  meshes,  in  this  chapter  we  investigate 
locking  on  curved  elements,  focusing  primarily  on  jo-refinement.  We  prove  that  on 
a  single  curved  element,  the  displacement  error  is  bounded  by  [k  —  a)~^  where  a 
depends  on  the  polynomial  degree  of  the  mapping  from  a  reference  element  to  the 
curvilinear  element.  We  adapt  the  mixed  methods  from  Chapters  3  and  4  to  curved 
elements  and  show  how  the  inf-sup  constant  can  be  estimated  when  these  methods 
are  used  on  curvilinear  meshes.  We  use  this  result  to  compute  inf  —  sup  constants  for 
several  mapped  elements  and  show  that  these  are  stable  for  practical  applications.  We 
also  show  that  modified  versions  of  MQl-6  can  be  proven  to  have  the  same  stability 
behavior  on  curvilinear  elements  as  on  parallelogram  elements. 

When  displacement  is  of  interest,  we  show  that  the  SFEM  shows  significant  de¬ 
terioration  when  curved  elements  are  used  on  the  boundary  with  p-refinement.  In 
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addition,  we  show  that  the  stress  error  is  adversely  affected  under  jo-refinement  by 
both  curved  boundary  elements  and  elements  with  curved  interior  edges.  We  also 
perform  experiments  to  show  that  the  mixed  methods  of  Chapters  3  and  4,  when  used 
on  a  curvilinear  mesh,  are  robust  for  a  model  problem. 

As  a  final  test,  in  Section  5.5,  we  use  a  benchmark  problem  from  [25]  to  investigate 
the  performance  of  SFEM  and  a  mixed  method  for  computing  point  values  of  stresses. 
In  [25],  the  SFEM  is  shown  to  be  ineffective  for  computing  point  values  of  the  sum  of 
normal  stress  and  a  post  processing  method  to  overcome  this  shortcoming  is  proposed. 
We  find  that  mixed  methods  with  curvilinear  meshes  can  be  used  to  accurately  extract 
these  point  values  without  any  post  processing. 

5.1  Single  Element  Analysis 

In  this  section  we  show  that,  on  a  single  element,  an  invertible  polynomial  mapping 
does  not  affect  the  asymptotic  rate  of  convergence  under  p-refinement,  but  can  affect 
the  relative  error,  resulting  in  a  locking  ratio  greater  than  one.  This  shift  is  shown  to 
depend  on  the  degree  of  the  polynomial  map.  This  result  is  given  for  the  case  that  K 
is  the  reference  square. 

Theorem  5.1  Let  Tk  :  K  K  be  an  invertible  polynomial  map,  J-k  €  [Qg{K)f 
and  Uk  e  [Vk{K)]'^,  where  Vk(K)  =  {v  :  v  =  v  o  ,  v  G  Qk(K)}.  Let  u  satisfy  (2.3) 


61 


on  K  and  U/^  he  the  solution  to  (2.4)-  If  K  is  the  reference  square  then 
(5-1)  ||u  —  UfcllHi(A')  <  C{k  — 

where  a  —  3q  —  2. 

Proof:  It  is  well  known  (see  Theorem  3.3  of  [1])  that  (5.1)  follows  if  we  can  show  it 
for  the  completely  incompressible  case.  Thus  we  consider  u  such  the  div  u  =  0,  and 
let  (j)  E  H^{K)  be  such  that  u  =  curl  solves  the  limiting  case  (A  — >■  oo)  of  (2.3). 

Let  the  globally  invertible  polynomial  map  J-fc  :  K  K  be  defined  as  y  =  Tk'^ 
where  Tk  €  [Qq{il)V'-  Let  Jk  he  the  Jacobian  of  J-k  and  d  be  the  determinant  of 
Jk-  Then 

^  ~  1^7 ^  ^2g-i,2g-2  U  'P2q-2,2q-i  C  Q2q-i-  ^  Since  iFx  IS  smooth  and 
invertible  on  if,  |d|  >  c  >  0  on  A'  and  for  <j)  G  H^{K)  we  have 

(5-2)  Ci\\(f)\\Hr^K)  <  W^WmiK)  ^  C2\\4>\\h-{K), 

where  ^  =  cj)  o  Tk-  Define  ^  For  any  I  there  exists  G  Qi{K)  such  that  for 

0<f<2,  r>l  (see  [21]) 

(5.3)  W'^  ~ *^\\'^\\h’'+^k)- 

Taking  and  combining  (5.2)  and  (5.3)  with  t  =  2,  we  have 

(5-4)  Wcf)  -  (f)k\\H^K)  <  C\\^  -  ^k\\m{K) 

^We  define  Pki,k2  =  spanfx^'a;^^  :  0  <  ui  <  ki,  0  <  a2  <  ^2}- 
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Now 


curl  (pk 
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Similarly,  every  other  term  in  curl  (pk  is  in  the  same  space.  Letting  a  =  Sq 
k  =  I  +  3q  —  then  I  =  k  —  a.  Then  (5.4)  implies 


(5.5) 


u  —  <  (^(fc  —  a)  ^Ilu||iyr(^p 


2  and 


□ 
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Theorem  5.1  gives  a  =  1  in  the  case  of  bilinear  mappings.  This  has  been  experiment¬ 
ally  observed  in  e.g.  Figure  5.2  of  [3],  where  the  shift  with  other  mappings  is  also 
computationally  investigated. 

Similar  results  are  available  if  K  is  the  reference  triangle.  In  that  case  if  J-k  € 
Vq{k)  and  Ujt  G  \yk{K)Y  where  Vk{K)  =  {u  :  u  =  u  o  G  Vk{k)},  then  a  = 

max{0, 3^  —  4}.  We  have  considered  the  simple  case  of  a  single  element  to  avoid  the 
complexities  of  more  general  meshes.  Experimentally,  the  single  element  results  carry 
over  qualitatively  to  the  multiple  element  case. 

5.2  Curvilinear  Mesh  Analysis 

We  now  consider  meshes  constructed  with  mapped  elements.  We  show  that  The¬ 
orem  3.1  applies  to  MQl-6  and  MTl-2  on  regular  (defined  below)  curvilinear  meshes 
provided  condition  (LS)  holds.  We  show  how  to  estimate  the  inf  —  sup  constant  (in 
terms  of  k)  for  these  mixed  methods  on  mapped  triangular  and  quadrilateral  elements 
for  certain  types  of  maps.  We  show  that  several  standard  mappings  produce  stable 
methods  and  we  provide  a  technique  for  examining  the  stability  of  additional  maps. 

Throughout  this  section,  we  assume  we  are  given  a  family  of  partitions  {Ch}  of  fl, 
where  fi  =  UKec^J^-  We  assume  the  mappings  Tk  '■  K  K  are  smooth  and  globally 
invertible  on  K.  We  let  Jk  denote  the  Jacobian  of  Tki  and  assume  its  inverse 
exists  for  any  x  £  K.  We  say  that  the  curvilinear  mesh  {Ch\  is  regular  [7]  if,  in 
addition  to  the  conditions  in  Section  3.1,  the  following  conditions  hold. 
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i)  Let  \  Jk\  =  IdetJft'l,  then 

(5.6)  ll’^'lloo,  ll-^A'^lloo,  sup|JAr(^)|,  [inf  |Ja(^)|]~^  <  a:, 

xeK 

where  /C  is  a  constant  independent  of  h. 

ii)  The  mesh  admits  the  existence  of  a  “Clement  operator”  [8];  a  continuous  inter¬ 
polate  using  averages  of  v  G  instead  of  point  values. 

Regular  meshes  can  be  constructed  on  many  domains  of  interest. 

To  use  Theorem  3.1  for  our  mixed  methods  on  curvilinear  meshes,  we  remark  that 
although  (Al)  and  (A2)  hold,  it  is  not  known  how  to  establish  (A3)  analytically.  We 
therefore  use  condition  (LS),  by  estimating  the  local  stability  condition  computation¬ 
ally.  In  addition  we  must  verify  that  (A4)  applies  for  these  mixed  methods.  The 
following  lemmas  establish  (A4)  for  the  case  of  triangular  and  quadrilateral  elements 
respectively. 

Lemma  5.1  Let  {Ch}  be  a  family  of  regular  ‘partitions  of  U,  and  let  Tk  '■  K  K 
denote  the  smooth  invertible  map  associated  ivith  K ,  with  K  the  reference  triangle. 
Then  the  pair  ,  'Po(fl))  is  globally  stable  where  Vk{Ll)  =  {u  ;  v\k  =  voTf^  .,v  G 

mk)}. 

Proof:  We  let  V  =  [JTo(n)]^,  Yh  =  ['P2{^)Y,  and  Qh  =  Vo{Ll).  We  use  Proposition 
4.1  from  [6],  which  says  it  is  sufficient  to  construct  an  operator  112  :  V  — >■  Yh  that  has 
the  following  two  properties: 

(5.7)  /  div  (v  —  Il2'v)qhdx  =  0  Vv  G  V  and  V  qn  &  Qh 

J 


and 
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(5.8)  ||n2  v||i,x  <  C{h],^Mo,K  +  |v|i,^)  Vv  G  V  and  €  C,. 

We  define  112  Vj/^  as  follows: 

U2^r^Ke[V2{K)]\ 

II2  V|/<'(M)  =  0  VM  =  vertex  of  K, 

(5.9)  II2  'V\K  ds  =  L  V  ds  Ve  =  edge  of  K. 

Using  (5.9),  we  establish  (5.7)  as  follows: 

div(v  -  112  v)  da:  =  -  ^2  'V\K)qh\K  dx 

KeCh 

=  XI  ^HK 

KeCh 

where  we  have  used  the  fact  that  qh  is  constant  over  each  K. 

Now  we  establish  (5.8).  First,  n2(v)  =  0  Vv  G  \Po[K)\^.  Also,  ||n2v||i  <  C'||v||i 
due  to  (5.9).  Hence  by  the  Bramble-Hilbert  lemma  we  have  |n2v]  <  C\v\i^k-,  or 
||n2v||o,ic  <  C\v\i^K-  Also,  by  an  inverse  estimate,  |n2v|i^/^  <  /i^^||n2v||o,A'-  Then 

||n2v||i,/c  <  ||n2v||o,7C  +  |n2v|i,/^ 

<  Civil, A- +  |n2v|i,A 

<  C|v|i,A  +  h^ln2V||o,A 

<  C'|v|i,A  +  ChX^^||v||o,A 

<  C'(/l^^l|v||o,A  +  |v|i,a). 


/  (v|A  -  n2  V|a)  •  nds  =  0, 
JdK 
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□ 

Lemma  5.2  Let  {Ch}  be  a  family  of  regular  'partitions  of  Cl,  and  let  Tk  :  K  ^  K 
denote  the  smooth  invertible  map  associated  with  K ,  with  K  the  reference  square. 
Then  the  pair  ([Q;c(0)]^,'Po(f^))  is  globally  stable  where  Qk{Cl)  =  {v  :  V|x  =  vo 

r^\-v€Q,(K)]. 

The  proof  is  similar  to  the  proof  of  Lemma  5.1. 

We  see  therefore  that  (A4)  holds  for  curvilinear  quadrilateral  and  triangular  regular 
elements.  Then  if  (LS)  holds,  the  inf-sup  constant  is  given  by  Theorem  3.1. 

We  investigate  condition  (LS)  computationally  as  follows.  Let  Tk  be  an  arbit¬ 
rary  smooth  invertible  mapping.  To  verify  the  local  stability  condition  (LS)  we 
should  use  the  subspaces  Yk{K)  H  Ho(/if)  and  Wk{K)  H  Ll{K).  However,  instead 
of  Wk{K)  n  Ll{K),  it  turns  out  to  be  easier,  both  computationally  and  theoretically, 
to  use  Wk{K)  n  Lq{K)  where 

Ll{K)  =  {V’  :  V’  -  ^  o  JT-i,  ^  e  Llik)}. 

Note  that  functions  ip  G  Lq{K)  satisfy  fj^  ip  |  dx  =  0  instead  of  ipdx  =  0.  We 
therefore  define  (LS),  an  alternative  condition  to  (LS): 

{LS)  For  every  K  ^  Ch  ,  given  q*  G  Wk{K)  D  Ll{K),  there  exists  v*  €  Yk{K)  fl 
[Hq{K)]‘^  such  that  (3.4)  holds. 

The  above  computational  substitution  can  be  justified  on  the  basis  of  the  following 


theorem. 
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Theorem  5.2  Let  the  spaces  Vjv,tTjv  satisfy  either  (LS)  or  (LS),  and  (A^)-  Then 
the  global  stability  estimate  (3.5)  holds. 

Proof:  The  case  (LS),  (A4)  is  proven  in  [20]. 

Let  q  e  Wn  be  arbitrary.  Then  for  each  K  ^  Ch,  q\K  —  qo^K  where  q  = 
q-\-  q*  ^  Wk{K),  with  q  G  VoiK)  and  fj^cfdx  =  0.  Then  q\K  =  (?  +  ?*)  o  = 
q  o  +90  =  q\ji  +  where  q\K  G  7^o(A^),  q*j^  G  Wk{K)  n  Ll{K).  Then  write 

q  =  q-\-  q*  with  q  being  the  weighted  projection  of  q  onto  the  space  of  piecewise 
constants: 

WN  =  {qeWN:  q\K  e  Po{K)  VA'  G  Ch}. 

From  Lemmas  5.1  and  5.2,  there  is  v  G  Vjv  such  that 

(5.10)  (divv,g)  >  CzWqWl  and  |v|i  <  C'4||g||o, 

with  positive  constants  C3  and  C4  independent  of  v  and  q. 

For  each  K  G  Ch,  we  have  q*^^  £  Wk{K)  D  Lq[K)  and  hence  by  {LS)  we  can  find 
V*  G  Vjv  such  that  G  [i7o(A")]^ 

(5.11)  (divv*,g*)  >  Cifc"'^||g*||o  with  |v*|i  <  C'2||g*||o- 
Since  v*^-  G  [Hq{K)Y,  it  holds  that 

(5.12)  (divv*,  ^)  =  0. 

Let  now  v  =  +  v*.  Using  (5.10)-(5.12),  the  Schwarz  inequality  and  the 
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arithmetic-geometric  mean  inequality  we  get 

(div  V,  q)  =  5(div  v,  q)  +  <5(div  v,  q*)  +  (div  v*,  q)  +  (div  v*,  q*) 

>  -  5|vMk*||o  + 

>  -  {SMUWl  -  {8el2)\v\l  +  C.k-'^hX 

>  8{C^  -  {C,em\q\\l  +  -  {8M)\\<f\\l 

>c.k-x\q\\l  +  hX) 

=  C,k-'^\\q\\l 

where  we  first  choose  e  =  C^IC^  and  then  8  =  ek~'^  fCi-  Since  we  then  also  have 
|v|i  <  <j|v|i  +  |v*|i  <  8C4\q\\o  +  C2\\q%  <  C\\q\\o, 
the  assertion  is  proved.  □ 

Since  (A4)  holds  by  Lemmas  5.1  and  5.2,  the  global  stability  condition  (3.5)  depends 
on  either  (LS)  or  (LS).  We  have  not  established  (LS)  for  curvilinear  elements  based 
on  the  mixed  methods  of  Chapters  3  and  4,  but  we  provide  computatioal  validation 
of  (LS)  in  the  next  section.  We  also  verify  (LS)  analytically  for  modified  versions  of 
our  quadrilateral  mixed  method  spaces  in  Section  5.4 

Provided  the  global  stability  condition  holds,  the  following  theorem  gives  error 
bounds  for  the  mixed  methods. 


Theorem  5.3  Let  Vjv(fi),  Wiv(n)  be  defined  on  regular  meshes  {Ch}  consisting  of 
curvilinear  triangles  and  quadrilaterals.  Let  the  global  stability  condition  (3.5)  hold 
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with  constant  7  independent  of  h.  Let  (A5)-(A6)  hold  with  n  >  0  in  (A6).  Then  for 
any  e  >  0,  there  exists  a  constant  Cg  independent  ofu,p,h,  and  k,  but  dependent  on 
tC  in  (5.6),  such  that 

(5.13)  ||u  -  unWi  <  +  \\p\\s-i) 

(5.14)  Up  -  pivllo  <  +  ||p||s-i) 

where 


I  =  min  {s  —  1,  A;  —  n}. 

See  proof  of  Theorem  5.2  in  [20].  For  MQl-4,6  and  MT2,  n  =  0.  For  MQ5  and  MTl, 

n  =  1. 

5.3  Computation  Of  Inf-Sup  Condition 

We  begin  with  the  results  of  computing  the  inf  —  sup  constant  from  (LS)  for  the  space 
MTl  on  mapped  triangular  elements.  We  consider  elements  that  have  exactly  one 
curved  edge.  The  first  map  is  a  parabola,  see  Figure  5.1,  ranging  from  the  extreme 
concave  edge  with  parallel  tangents  to  a  reasonable  range  of  convex  parabolic  edges. 
The  results  are  found  in  Figure  5.2.  We  see  that  as  long  as  the  case  where  the  edges 
meet  with  parallel  tangents  is  avoided,  the  inf  —  sup  constant  has  about  the  same 
decay  rate  as  the  unmapped  reference  triangle,  where  7  RS  1 . 


Figure  5.2:  Inf-Sup  Constants  -  Parabolic  Triangles  -  Concave,  Convex 


The  next  case  is  for  the  circular  maps  shown  in  Figure  5.3.  The  range  of  maps 
considered  is  again  from  the  case  of  concave  parallel  tangents  to  a  reasonable  range  of 
convex  circular  edges.  As  with  parabolas,  as  long  as  the  case  where  edges  meet  with 
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Figure  5.3:  Curvilinear  Triangles  -  Circular  Edges 


Figure  5.4:  Inf-Sup  Constants  -  Circular  Triangles  -  Concave,  Convex 


parallel  tangents  is  avoided,  the  decay  rate  7  is  comparable  to  that  of  the  immapped 
element.  This  is  consistent  with  our  definition  of  regular  mesh,  which  requires  element 
interior  angles  be  bounded  away  from  0  and  tt.  These  results  are  seen  in  Figure  5.4. 


InfSup  Constant 
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Figure  5.5:  Curvilinear  Quadrilaterals  -  Parabolic  Edges 


Figure  5.6:  Inf-Sup  Constants  -  Parabolic  Quadrilaterals  -  Concave,  Convex 

We  conclude  that  if  triangular  meshes  are  constructed  so  as  to  produce  sufficiently 
regular  elements,  the  resulting  inf —  sup  constants  are  comparable  to  the  linear  case. 
We  note  that  automatic  mesh  generators  exist  that  produce  regular  meshes  for  most 
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Figure  5.7:  Curvilinear  Quadrilaterals  -  Circular  Edges 


Figure  5.8:  Inf-Sup  Constants  -  Circular  Quadrilaterals  -  Concave,  Convex 


We  performed  similar  calculations  for  the  curvilinear  quadrilaterals  in  Figures  5.5 
and  5.7,  using  MQ4.  The  results  are  given  in  Figures  5.6  and  5.8.  We  see  that  the  k 
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tested  are  well  within  the  pre-asymptotic  range,  since  the  rate  is  not  close  to  — 


5.4  A  Stable  Curvilinear  Quadrilateral  Element 


The  stability  of  mixed  methods  MQl-6  depends  on  condition  (LS)  (or  equivalently 
(LS))  and  the  inf-sup  constant  for  these  methods  depends  on  7  in  (LS).  So  far  we  only 
have  computational  estimates  of  7.  In  this  section  we  define  curvilinear  elements  for 
which  a  theoretical  inf-sup  constant  of  Ck~^  can  be  established,  identical  to  that  for 
parallelogram  elements  (see  Corollary  3.1). 

To  define  these  elements,  use  (3.2)  to  define  Vn,Wn  in  terms  of  Vk{K),Wk{K). 
For  MQl-6,  we  define  Wk{K)  by  (3.1).  We  define  Vk{K)  as  follows: 


(5.15) 


n(A')  = 


V 


^k(v) 


:  V  e 


V^(A0 


> 


9 


:veV,(^)\V0(A0  J 

where  u  =  =  Ua  I-OkU  e  is  the  Piola  transformation  of  u  [6].  Since 

V  e  v,\K)  ^  V  e  v«(K),  Vjv  is  again  conforming.  We  then  have  the  following 
theorem. 


Theorem  5.4  Let  V]v,IT}v  be  defined  by  MQl-6,  (5.15)  and  (3.2).  Then 
(5.16) 


(divv,q)  i 

inf  sup  I,  >Ck  2, 


where  the  constant  C  is  independent  of  h,k  and  q,  i.e.  (3.5)  holds. 


Proof:  We  show  that  (LS)  holds  with  7  =  |.  The  theorem  then  follows  from  The¬ 


orem  5.2. 
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Let  q*  G  Wk{K)  H  Lq(K)  be  arbitrary  and  define  q*  G  Wk{K)  H  Lq(K)  by  q*  = 
q*oFK  (q*  ^  ^oi^)  definition  of  Lq{K)  ).  On  the  reference  square,  the  continuous 
inf-sup  condition  holds.  Hence,  there  exists  v  G  [Hq{K)Y  such  that 

(5.17)  (divv,q*)  >  (7119*11^,  and  <  ||q*||o_^, 

with  C  independent  of  q*  and  v. 

Define  v*  G  V°(^)  by  v*  =  T^v.  By  (A2),  (A3)  and  (5.17)  we  obtain  (integrating 
by  parts) 

(5.18)  (div  V*,  q*)j^  =  -(v*,  vr)^  =  -(TfcV,  S7Q*)k 

=  -(v*,Vfk  =  (divv,q*)^ 

>  cwrwiK 

and 

(5.19)  =  ir.vl,,^  <  CA:^||v||,,^.  <  Ck^r\\o,K- 

This  proves  the  discrete  inf-sup  condition  on  the  reference  square. 

Next,  we  let 

(5.20)  V*  =  CK{r)  =  W-^JK-r  O  Tk\ 

Then  by  (Al),  £k(v*)  G  Yk{K)  fl  [Hq(K)Y.  Using  the  basic  property  of  the  Piola 
transformation  [6],  we  have 

(div  V*,  q*)K  =  (div  v*,  q*)^. 


(5.21) 
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By  assumption  of  the  reqularity  of  the  elements,  we  have 

(5.22)  Ik’lloA  SS  c4lirilo,;f  ““d  |v'li,K  <  C'fc7|v*l,,K. 

Then  {LS)  follows  from  (5.21)  and  (5.22),  and  the  local  conditions  (5.18)  and  (5.19). 

□ 

By  Theorem  5.3,  since  (3.5)  holds  and  (A5)-(A6)  hold,  therefore  we  get  the  conver¬ 
gence  rates  in  (5.13)  and  (5.14).  We  have  not  used  these  modified  spaces  in  our 
computations.  However,  the  experiments  in  Section  5.5  reflect  these  rates  of  conver¬ 
gence. 

We  point  out  that  to  implement  these  spaces,  one  needs  only  a  basis  for  Wk{K) 
and  Vk{K),  the  mapping  Tk  and  its  Jacobian.  In  particular,  the  inverse  mapping 
is  not  required. 

5.5  Numerical  Experiments 

In  the  p-version,  element  size  is  fixed,  so  curvilinear  elements  must  be  used  to  accur¬ 
ately  model  the  domain  for  most  practical  problems.  We  use  the  blending  function 
method  (see  [24])  to  match  the  domain  exactly.  We  also  study  meshes  which  have 
interior  edges  that  are  curved.  We  study  the  standard  and  mixed  methods  on  both 
curvilinear  quadrilateral  and  triangular  elements. 

We  limit  our  investigation  to  domains  (meshes)  that  are  typical  in  applications. 
This  is  of  course  a  subjective  restriction,  but  the  curved  domains  seen  in  practice 
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QDIST 


1 


QPARA  QCIRC 


Figure  5.9:  Curvilinear  Quadrilateral  Meshes 


can  often  be  parameterized  using  piecewise  low  degree  polynomials  or  trigonometric 
functions  (ellipses).  These  domains  can  be  produced  by  mapping  the  reference  domain 
with  non-affine  invertible  mappings. 

The  quadrilateral  meshes  we  use  are  given  in  Figure  5.9  and  the  triangular  meshes 
are  given  in  Figure  5.10.  These  meshes  were  chosen  to  include  elements  with  edges  that 
are  parameterized  by  low-order  polynomials  (parabolas)  and  trigonometric  functions 
(circles),  with  both  interior  and  boundary  curved  edges. 

5.5.1  Finite  Element  Spaces 

In  Chapters  3  and  4  we  described  several  polynomial  spaces  defined  on  reference 
elements.  When  curvilinear  elements  are  used,  the  finite  element  spaces  are  com¬ 
posed  by  applying  non-affine  invertible  smooth  mappings  to  reference  polynomial 
spaces.  The  resulting  spaces  are  not  polynomials,  but  mapped  images  of  these  ref¬ 
erence  polynomial  spaces.  Recall  that  SQ  :  14  =  Q'ki^)  ST:  14  =  Vk{K). 


Figure  5.10:  Curvilinear  Triangular  Meshes 


We  restrict  our  attention  to  MQl:  14  =  Qk{K)  H  VkJr2{K)^Wk  =  Vk-i{K)  and 
MTl  :  14  =  Vk{K)i  Wk  =  Vk-2{K).  A  partition  Ch  of  the  domain  0  is  constructed 
such  that  fl  =  UxeChK  as  indicated  in  Figures  5.9  and  5.10.  For  each  K  there  exists  a 
smooth  invertible  mapping  such  that  K  =  Tk{K),  where  K  is  either  the  reference 
triangle  or  square. 


5.5.2  Computational  Convergence  Results 

We  restrict  our  attention  to  the  p  version  and  begin  by  examining  SQ  and  ST  on 
meshes  with  curvilinear  quadrilateral  and  the  triangular  elements,  respectively.  We 
continue  using  energy  norm  and  SNS  to  measure  the  error  as  discussed  in  Chapter 
3.  Throughout  this  section  we  solve  model  problem  SING  with  a  near-singular  (un¬ 
smooth)  solution. 
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Degrees  of  Freedom  —  k  =  2 . 10  Degrees  of  Freedom  —  k  =  2,  10 

Figure  5.11:  Method  SQ,  Square  Mesh  vs  QDIST,  Unsmooth  solution. 

SFEM  on  Quadrilateral  Elements 

We  first  look  at  method  SQ  on  meshes  with  straight  and  curved  interior  boundaries. 
For  this  we  solve  problem  SING  on  meshes  QDIST  and  QPARA  and  also  a  mesh  with 
four  square  elements.  In  Figures  5.11  and  5.12  we  see  that  the  error  in  energy  norm  is 
not  sensitive  to  changes  in  i/,  however  this  is  not  true  for  the  SNS  error.  As  i/  — |  the 
relative  error  in  SNS  increases  significantly  on  both  QDIST  and  QPARA,  much  more 
so  than  it  does  for  the  square  mesh.  We  conclude  that  method  SQ  has  near  optimal 
behavior  for  energy  norm  but  that  it  locks  in  SNS  when  the  mesh  contains  elements 
with  curved  interior  edges.  Table  5.5.2  shows  that  the  locking  ratios  are  significantly 
worse  when  curved  interior  edges  are  present. 

The  results  for  the  model  problem  on  mesh  QCIRC  are  shown  in  Figure  5.13.  We 
see  that  the  error  in  both  the  energy  norm  and  SNS  increases  significantly  as  z/  — 
Thus  we  conclude  that  the  standard  method  degrades  significantly  when  curvilinear 
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N 

95 

135 

183 

239 

303 

375 

455 

543 

Square  -  Ls 

7.12 

8.85 

7.85 

6.71 

5.63 

4.33 

3.19 

2.53 

QDIST  -  Ls 

12.08 

45.26 

28.78 

57.90 

61.98 

66.20 

64.77 

65.72 

QPARA  -  Ls 

35.93 

23.70 

57.73 

64.52 

65.49 

74.49 

78.79 

81.54 

Table  5.1:  SNS  Locking  Ratios  (L5(0.4999,  A^)),  SQ,  Unsmooth  Solution. 


Figure  5.12:  Method  SQ,  Square  Mesh  vs  QPARA,  Unsmooth  solution. 


elements  are  used  to  partition  domains  with  curved  boundaries.  This  energy  norm 
locking  was  not  observed  for  method  SQ  on  straight  sided  elements. 


SFEM  On  Triangular  Elements 

To  evaluate  the  standard  method  on  curvilinear  triangular  elements  we  use  meshes 
TPARA  and  TCIRC  shown  in  Figure  5.10.  Mesh  TPARA  has  curved  interior  edges, 
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Figure  5.13:  Method  SQ,  Mesh  QCIRC,  Unsmooth  solution. 


Figure  5.14:  Method  ST,  Mesh  TPARA,  Unsmooth  solution. 


while  TCIRC  has  curved  elements  on  the  boundary.  As  seen  in  Figures  5.14  and  5.15, 
method  ST  exhibits  no  energy  norm  locking  on  either  mesh  while  SNS  locking  ratios 
are  significant  in  both,  as  seen  in  Table  5.5.2. 


nu  «  0.4999 
nu  =  0.3 


^  -  nu  »  0.4999 
o  -  nu  =  0.3 


Degrees  of  Freedom  —  k  =  2  ...  1 1 


10  10 

Degrees  of  Freedom  —  k  =  2  ... 


Figure  5.15;  Method  ST,  Mesh  TCIRC,  Unsmooth  solution. 


47 

95 

159  1 

1 

239 

335 

447 

575 

719 

879 

20.82 

26.10 

30.50 

34.10 

37.40 

40.72 

44.04 

47.32 

50.29 

N 

35 

71 

119 

179 

251 

335 

431 

539 

659 

39.55 

37.26 

32.07 

15.73 

61.06 

136.21 

84.67 

46.92 

123.77 

Table  5.2:  SNS  Locking  Ratios  (L5(0.4999,  A^)),  ST,  Unsmooth  Solution. 


MFEM  On  Quadrilateral  Elements 

We  now  turn  our  attention  to  mixed  finite  element  methods.  We  do  the  same  sequence 


of  experiments  for  the  mixed  methods  as  we  did  for  the  standard  method.  When  curved 
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Figure  5.16:  Method.  MQl,  Mesh  QDIST,  Unsmooth  solution. 


Figure  5.17:  Method  MQl,  Mesh  QPARA,  Unsmooth  solution. 


edges  are  in  the  mesh  interior,  as  in  QDIST  and  QPARA,  the  curves  for  u  =  0.3  and 
for  V  =  0.4999  are  essentially  the  same.  Thus  we  see  in  Figures  5.16  and  5.17  that 
method  MQl  is  very  robust  on  meshes  with  curved  interior  edges.  Similarly,  on  mesh 
QCIRC  virtually  no  increase  in  the  errors  is  observed,  as  seen  in  Figure  5.18. 
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Figure  5.18:  Method  MQl,  Mesh  QCIRC,  Unsmooth  solution. 


Figure  5.19:  Method  MTl,  Mesh  TPARA,  Unsmooth  solution. 


MFEM  On  Triangular  Elements 

To  study  mixed  methods  on  curvilinear  triangular  elements  we  use  MTl  (which  has 
the  same  displacement  space  as  ST)  to  solve  the  model  problem  on  mesh  TPARA  as 
an  example  of  a  mesh  with  curved  interior  edges  and  on  mesh  TCIRC  as  and  example 
of  a  mesh  with  curved  edges  on  the  boundary.  We  see  in  Figures  5.19  and  5.20  the 
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Figure  5.20:  Method  MTl,  Mesh  TCIRC,  Unsmooth  solution. 

error  is  unaffected  as  u 

In  summary,  when  the  curved  edges  are  on  the  mesh  interior,  the  standard  method 
showed  mild  energy  norm  locking  but  significant  SNS  locking  on  both  quadrilat¬ 
eral  and  triangular  meshes.  On  the  other  hand,  the  mixed  methods  show  essentially 
identical  error  curves  for  both  quadrilateral  and  triangular  meshes  in  both  energy  and 
SNS  as 

Meshes  with  curved  exterior  boundaries  influence  both  displacement  and  SNS 
locking,  for  the  standard  method.  Some  locking  was  seen  in  the  energy  norm  on 
both  quadrilateral  and  triangular  meshes  which  was  not  present  on  straight  sided 
elements.  In  addition,  SNS  locking  was  severe.  For  the  mixed  methods  the  rate  of 
convergence  was  unaffected  in  the  problems  we  studied,  both  for  the  energy  norm  and 
SNS,  although  minor  locking  ratios  were  observed  in  some  cases. 
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5.6  Computing  Point  Values  -  A  Benchmark  Prob¬ 
lem 


In  the  engineering  design  process,  often  the  objective  is  to  identify  regions  of  high 
stress  and  produce  a  design  that  distributes  the  stress  more  evenly  throughout  the 
body,  producing  lower  maximum  stress  levels.  This  requires  the  analytical  tool  to 
accurately  compute  point  values  of  stress  and  other  quantities  of  engineering  interest. 
To  evaluate  the  point  value  extraction  capability  of  the  FEMs  we  are  studying,  we 
turn  to  a  benchmark  problem  from  [25]  which  we  will  refer  to  as  the  Rigid  Circular 
Inclusion  (RCI)  problem. 


Figure  5.21:  Rigid  Circular  Inclusion  Problem  and  Quadrilateral  Mesh. 

This  problem  involves  a  rigid  circular  inclusion  in  an  infinite  plate,  subject  to 
unidirectional  tension  under  plane  strain  conditions.  The  problem  is  illustrated  in 
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Figure  5.22:  Rigid  Circular  Inclusion  Problem  -  Triangular  Mesh. 
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Figure  5.23;  RCI  problem,  Method  SQ,  Quadrilateral  elements. 


Figure  5.24:  RCI  problem,  Method  ST,  Triangular  elements. 

Before  we  look  at  SNS  point  value  extraction,  we  do  the  same  error  measurements 
we  have  done  in  previous  chapters.  The  first  two  experiments  are  with  SQ  and  ST 
using  p-refinement  on  the  RCI  problem;  the  triangular  mesh  is  shown  in  Figure  5.22. 
The  results  are  shown  in  Figures  5.23  and  5.24.  In  Table  5.6  we  see  that  both  methods 
have  significant  locking  ratios,  although  the  rates  of  convergence  appear  unaffected  as 


%  Rel  Error  -  Energy  Norm 
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z/  In  general,  the  locking  ratios  for  ST  are  smaller  than  those  for  SQ.  In  both 

methods,  the  locking  ratios  stay  relatively  constant  over  the  range  of  k  we  tested. 


SQ 

N 

40 

69 

96 

136 

184 

240 

304 

376 

Le 

11.21 

7.52 

9.14 

8.84 

8.67 

9.95 

15.05 

32.24 

Ls 

585.28 

296.40 

326.04 

303.42 

318.91 

329.10 

382.76 

377.50 

ST 

N 

B 

54 

96 

150 

216 

294 

384 

486 

Le 

19 

3.60 

3.59 

3.75 

3.63 

3.62 

3.67 

3.46 

Ls 

136.61 

120.29 

183.64 

178.30 

93.06 

82.93 

129.70 

159.59 

Table  5.3:  Locking  Ratios,  SQ  and  ST,  RCI  Problem. 


Degrees  of  Freedom  Degrees  of  Freedom 


Figure  5.25:  RCI  problem.  Method  MQl,  Quadrilateral  elements. 


90 


Degrees  of  Freedom  —  k  =  2  ...  10  Degrees  of  Freedom  —  k  =  2  ...  10 


Figure  5.26;  RCI  problem,  Method  MTl,  Triangular  elements. 

In  both  energy  norm  and  SNS,  MQl  and  MTl  produce  locking  ratios  of  about  2, 
as  seen  in  Figures  5.25  and  5.26. 

We  conclude  that  mixed  methods  on  meshes  with  curvilinear  boundaries  are  sub¬ 
ject  to  slight  increase  in  error  as  i/  ^  although  the  rate  of  convergence  seems 
unaffected.  On  the  other-hand,  SFEM  on  meshes  with  curvilinear  boundaries  shows 
significant  locking  ratios,  particularly  in  SNS.  For  this  model  problem,  we  also  see 
non-trivial  locking  ratios  in  energy  norm,  which  are  significantly  greater  than  those 
observed  for  the  model  problem  of  Section  5.3. 

As  mentioned  above,  an  important  capability  of  an  FEM  is  point  extraction  of 
stress,  which  in  SQ  and  ST  may  be  done  by  differentiating  the  computed  displace¬ 
ments.  It  is  shown  in  [25]  that  —  <Ty  and  r^y  can  be  accurately  extracted  this  way, 
even  as  ->  |,  but  that  ax  +  (Ty  —  SNS  cannot.  (For  the  latter,  at  least  in  two 
dimensions,  post-processing  could  be  used  -  see  [25].) 
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Figure  5.27:  SNS  Point  Values  MQl  vs  SQ,  k  =  3,  =  0.3  and  0.4999. 


Figure  5.28:  SNS  Point  Values  MQl  vs  SQ,  k  =  8,  z/  =  0.3  and  0.4999. 

In  Figure  5.27  and  5.28  we  compare  the  point  values  of  SNS  along  the  arc  AE  (see 
Figure  5.21),  computed  using  SQ  (differentiation  of  the  displacements)  and  MQl.  We 
note  that  when  v  —  0.3,  both  methods  compute  the  actual  SNS  reasonably  well  even 
for  k  =  3]  while  for  fc  =  8  they  match  the  actual  values  very  well.  On  the  other  hand, 
when  1/  =  0.4999,  we  observe  that  even  for  =  3,  MQ3  still  matches  the  exact  solution 
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Figure  5.29:  SNS  Point  Values  MTl  vs  ST,  k  =  3,  u  =  0.3  and  0.4999. 


Figure  5.30:  SNS  Point  Values  MTl  vs  ST,  k  =  8,  u  =  0.3  and  0.4999. 


reasonably  well,  while  SQ  produces  highly  oscillatory  results.  The  oscillations  continue 
in  SQ  even  for  k  =  8,  while  MQ3  and  the  exact  solution  are  essentially  identical  for 
this  k. 

We  performed  the  same  set  of  experiments  with  ST  and  MTl  on  a  triangular  mesh 
(see  Figure  5.21)  for  RCIP  with  6  elements.  The  results  in  Figures  5.29  and  5.30  are 
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qualitatively  the  same  as  with  quadrilateral  elements.  MTl  performs  well  over  the 
range  of  v  tested  while  ST  produces  good  results  for  v  =  0.3  but  when  v  =  0.4999 
the  output  from  ST  is  oscillatory  even  for  large  k. 

5.7  Summary 

In  this  chapter,  we  have  investigated  the  locking  characteristics  of  several  different 
FEMs.  The  performance  of  these  methods  is  summarized  in  Table  5.7.  A  “Yes” 
means  that  the  method  converges  at  the  asymptotic  best  rate,  is  locking-free,  and  has 
negligible  locking  ratios  (<  4).  “Rate”  means  the  rate  of  convergence  is  less  than 
optimal.  “Lock”  means  the  method  locks  to  some  degree.  “LR”  means  the  locking 
ratios  are  significant  (>  4).  The  h  version  results  refer  to  experiments  on  meshes  with 
straight  edges  (Chapter  3  and  4  results).  The  p  version  results  refer  to  meshes  with 
straight  edges  as  well  as  meshes  with  curved  boundaries  (Chapter  5  results). 

No  computations  were  done  using  MT2,  although  we  expect  it  to  perform  similar 


to  MQ1,2  and  4. 


Method 

h  version 

p  version 

Displacement 

Stress 

Displacement 

Stress 

Straight 

Straight 

Straight 

Curve 

Straight 

Curve 

SQ 

Lock 

Lock 

Yes 

Lock 

LR 

Lock 

ST 

Lock  fc  <  3 

Lock  A;  <  3 

Yes 

LR 

Yes 

LR 

Yes,  k  =  4 

LR  A;  =  4 

MQ1,2,4 

Yes 

Yes 

Yes 

Yes 

Yes 

Yes 

MQ3,5 

Rate 

Rate 

Yes 

Yes 

Yes 

Yes 

MTl 

Rate 

Rate 

Yes 

Yes 

Yes 

Yes 

Table  5.4:  Convergence  and  locking  properties  of  the  different  methods. 


Chapter  6 


An  Application  to  Nonlinear 
Elasticity 


The  analysis  of  mechanical  components  and  structures  is  most  commonly  performed 
on  the  basis  of  the  linear  theory  of  elasticity  [12].  This  approach  gives  excellent  results 
for  a  large  variety  of  problems  and  is  based  on  the  assiunption  that  the  strains  and 
displacements  are  small.  In  the  case  of  strains,  “small”  means  much  less  than  unity.  In 
the  case  of  displacements,  “small”  means  that  the  original  undeformed  configuration 
and  the  final  configuration  of  a  body  should  not  differ  by  much,  so  that  enforcing 
equilibrimn  in  the  imdeformed  configuration  is  the  same  as  enforcing  equilibritim  in 
the  deformed  configuration.  In  many  important  problems,  for  example  when  the 
limits  of  the  load  carrying  capacity  of  a  structure  are  of  interest,  the  displacements 
are  large  enough  that  the  equilibrium  equations  in  the  deformed  configuration  are 
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very  different  from  those  in  the  undeformed  configuration.  In  order  to  obtain  reliable 
results,  a  geometrically  nonlinear  analysis  must  be  performed  [16]. 

In  [16],  an  iterative  method  is  described  for  solving  geometrically  nonlinear  prob¬ 
lems  in  finite  elasticity.  It  is  based  on  a  spatial  formulation,  which  has  some  advantages 
over  the  traditional  material  formulation.  These  include:  (a)  the  equilibrium  equa¬ 
tions  are  satisfied  in  the  deformed  configuration;  (b)  the  displacement  components  are 
approximated  by  high-order  polynomials  and  the  blending  function  method  is  used  to 
provide  an  accurate  geometric  description  of  the  deformed  configuration;  (c)  the  dis¬ 
cretization  errors  are  controlled  by  p-refinement;  (d)  the  simulation  of  nonconservative 
loads,  such  as  follower  loads,  does  not  require  any  additional  extensions  in  the  formu¬ 
lation,  nor  does  it  destroy  the  symmetry  of  the  resulting  stiffness  matrix;  and  (e)  the 
extraction  of  stresses  from  the  finite  element  solution  is  straightforward,  since  in  this 
approach  the  equilibrium  equations  are  expressed  in  terms  of  the  Cauchy  stresses. 
This  formulation  is  a  standard  finite  element  displacement-based  formulation.  We 
use  this  spatial  formulation  to  define  a  mixed  method  by  introducing  an  independent 
variable  for  the  “pressure”. 


6.1  Mixed  Method  -  Spatial  Formulation 
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Using  the  equilibrium  equations  of  elasticity  (2.1a)  we  multiply  by  a  test  function 
V  €  V  and  integrate  to  get  the  variational  form  of  the  equilibrium  equations, 

(6.1)  /  (cr^Ui,i  +  (TyV2^2  +  +  U2,2))  dx  = 

J 

/  fiVi  dx+  giVi  ds  Vv  G  V. 

Ja  Jr, 


In  the  case  of  plane  strain 


2p  +  A  A  0 
A  2//  "b  A  0 


0  /i 


The  nonlinear  strain  tensor  e  is 


=  a  —  b 


where 


is  the  linear  part  of  e  and 


«1,2  +  ^^2.1 


|«i  +  «l.i) 

|(u2  2  +  u2  2) 


^1,1^1, 2  +  ^^2, 1^*2, 2 
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is  the  nonlinear  part.  Substituting  (6.2)  into  (6.1)  we  get 

(6.6)  /  (2/i(£a;Ui,i  +  £yU2,2)  +  M7x3/(t^i,2  +  V2,i)  +  +  ey)(i^i,i  +  ■^2,2))  dx 

J  O 

/  fiVi  dx+  QiVi  ds  Vv  G  V. 


Define 


(6.7)  p  = -A(£^ +£j,), 

and  substitute  (6.3)  and  (6.7)  into  (6.6).  We  move  the  nonlinear  terms  to  the  right 
hand  side  to  get 

(6.8)  2/2  /  (aiui,i  +  02^2,2  +  i  a3(ui.2  +  t^2,i))  dx  -  I  p(ui,i  +  ^2,2)  dx  = 

Ja  2  Jn 

jf  fiVi  dx  +  QiVi  ds  +  2/2  ^  (hiui.i  +  62^2,2  +  ^&3(^^i,2  +  ^^2,0)  dx. 

Putting  (6.7)  into  variational  form,  we  get 

(6.9)  /  (ui  +  a2)w  dx  +  X~^  I  pw  dx  = 

Jn  Ja 

f  [bi  +  62)22;  dx  V22;  G  W. 

J 

The  system  of  equations  (6.8)  and  (6.9)  is  the  variational  form  of  our  mixed  method. 
We  produce  an  iterative  scheme  as  follows.  Let  designate  the  reference  configur¬ 
ation.  Let  upper  indices  “(?)”  represent  quantities  associated  with  the  iteration. 
Then  we  have:  uh+^)  is  the  unknown  displacement  vector  to  be  determined  at  the 
iteration;  u^)  is  the  known  displacement  at  the  iteration;  and  Oh)  and  define  the 
domain  and  the  traction  boundary  of  the  deformed  configuration  at  the  iteration. 
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Using  (6.8)  and  (6.9)  we  define  the  iterative  procedure 

(6.10)  2fj,  +  a^2^^^V2,2  +  Vi,2  +  ^^2,i)) 

-  L  + ”’•=)  = L  + /<■) 

2/J  (M’*”!,!  +  4'*»2,2  +  j4’*('’1.2  +  ”2,1)) 


+ 


and 


(6.11)  /  +  ^2  J^^)^  ^  = 

jn(')  ’  ’  9a(’) 

^  X(0  ^ 

We  use  the  mixed  method  spaces  MQl  and  first  solve  the  discrete  version  of  (6.10) 
and  (6.11)  for  by  integrating  on  the  reference  configuration  which  is  the 
linear  problem  (uj^^  =  0).  This  produces  the  displacement  vector  such  that 
Using  we  again  solve  the  discrete  form  of  (6.10)  and  (6.11)  for 
by  integrating  on  Continuing  in  this  way,  we  produce  a  sequence  of  solutions 
{ujy  }  which  converges  to  u^,  the  finite  element  approximation  to  the  solution  to  (6.1). 
We  also  get  an  approximation  to  p. 

At  each  iteration  we  have  the  following  system  of  equations 


(6.12) 

(6.13) 


Au  —  Bp  =  Ri 


B^u  +  Cp  =  R2. 


On  each  element,  we  condense  out  the  pressure  by  solving  (6.13)  for  p  =  C  ^(R2  — 
B^u)  and  substituting  into  (6.12),  giving 


(6.14) 


(A  +  BC-^B^)u  =  Ri  +  BC-^R2. 
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Figure  6.1:  Cantilever  Subject  to  Pressure  Load 
We  then  solve  6.14  for  u.  By  saving  C~\R2  and  B  on  each  element,  we  can  recover 

P- 

Since  each  iteration  involves  integrating  over  a  mapped  domain,  each  element  in 
the  reference  mesh  becomes  a  curvilinear  element  with  potentially  no  straight  edges. 


6.2  Numerical  Results 


We  test  the  mixed  method  on  the  problem  illustrated  in  Figure  6.1.  A  beam  of  length 
I  and  height  h  is  fixed  on  one  end  and  subject  to  a  uniform  pressure  load  —  L/2  on 
top  and  L/2  on  bottom.  We  set  body  forces  f  =  0.  In  [16]  this  problem  and  several 
others  are  solved  using  the  SFEM  spatial  formulation.  Here,  we  are  interested  in  the 
case  where  u  ^  which  was  not  investigated  in  [16].  We  want  to  approximate  both 


the  displacements  and  the  sxun  of  normal  stresses,  SNS  which  is  related  to  p  as  before. 


— 2(/i  +  A) 


We  compute  the  displacement  and  SNS  for  the  case  when  /  —  10,  h  =  1,  modulus 
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of  elasticity  E  —  12000,  Poisson  ratio  u  =  0.3  and  0.4999  and  the  load  L  =  5.  We  use 
SFEM  with  Vk{K)  =  Qk{K)  and  our  mixed  method  with  MQl;  in  both  cases  we  use 
one  quadrilateral  element.  The  displacement  results  are  given  in  Table  6.1. 


Standard  Method 

Mixed  Method 

//  =  0.3 

u  =  0.4999 

//  =  0.3 

z/  =  0.4999 

-2.4590 

-0.2838 

-2.6317 

-2.4947 

-4.0682 

-1.1477 

-4.2534 

-3.6171 

fc  =  4 

-4.8943 

-2.0605 

-4.9601 

-3.9954 

k  =  5 

-5.0443 

-3.1704 

-5.0731 

-4.1209 

fc  =  6 

-5.0762 

-3.6285 

-5.0895 

-4.1787 

k  =  7 

-5.0949 

-3.8592 

-5.1037 

-4.2243 

k  =  8 

-5.1042 

-4.0128 

-5.1390 

-4.3200 

Table  6.1  Vertical  Displacements  of  Point  A. 


The  displacement  for  =  0.3  is  about  the  same  for  both  methods,  with  the  biggest 
difference  when  k  is  small.  However,  when  i/  =  0.4999,  the  difference  is  significant, 
especially  for  small  k.  Assuming  the  results  for  A;  =  8  are  reasonably  close  to  the 
true  solution,  we  see  that  the  mixed  method  produces  more  accurate  displacement 
approximations  when  k  is  small. 

We  are  primarily  interested  in  computing  point  values  of  SNS  when  the  material 
is  nearly  incompressible.  We  use  the  standard  formulation  to  extract  SNS  for  u  =0.3 
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and  0.4999.  The  results  are  given  in  Figure  6.2  for  various  values  of  k.  These  results 
are  qualitatively  the  same  as  those  in  Chapter  3,  where  the  standard  method  produced 
highly  oscillatory  SNS  point  values  as  i/  — We  note  that  the  scale  for  SNS  when  i/ 
—  0.4999  is  0(10^),  10  times  the  scale  when  i/  =  0.3.  On  the  other  hand,  the  mixed 
method  produced  the  results  shown  in  Figure  6.3.  Even  though  the  exact  solution 
is  not  known,  we  see  very  intuitive  distribution  of  SNS  along  the  beam  (with  no 
qualitative  change  as  i/  — >■  |);  with  a  maximum  at  x  =  0  smoothly  decaying  to  0  at 
the  end  (x  =  10)  of  the  beam. 


Figure  6.2;  SNS  on  Top  Edge,  Standard  Method 

We  conclude  that  mixed  methods  can  effectively  eliminate  stress  locking  in  some 


geometrically  non-linear  problems. 


£ 


Figure  6.3:  SNS  on  Top  Edge,  Mixed  Method 
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